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PREFACE 


A theoretical analysis solves for the steady-state photocur- 
rents produced by a given photogeneration rate function with 
negligible recombination in simple silicon diodes, consisting of 
a uniformly doped quasi-neutral region (called "substrate" below) 
between a p-n junction depletion region (DR) and an ohmic contact 
(electrode) . Special attention is given to conditions that pro- 
duce "funneling" (a term used by the single-event-effects commu- 
nity) under steady-state conditions. Funneling occurs when carri- 
ers are generated so fast that the DR becomes flooded and par- 
tially or completely collapses. Some or nearly all of the applied 
voltage plus built-in potential normally across the DR is now 
across the substrate. This substrate voltage drop affects sub- 
strate currents. The steady-state problem can provide some quali- 
tative insights into the more difficult transient problem. Chap- 
ter 6 discusses some similarities between the steady-state and 
transient problems. 

The DR boundary (DRB) is defined by an equation, but can be 
recognized from computer simulation results by plotting electron 
and hole densities, against a spatial coordinate, together on the 
same graph. Such a plot shows a reasonably well defined boundary 
that separates a space-charge region from a quasi-neutral region. 
With the DRB reasonably well defined, DR and substrate voltage 
drops are also reasonably well defined, and quantify the extent 
of DR collapse and the strength of funneling. A collapsed DR can 
also be recognized by a small width. 

It was found that the substrate can divide into two subregions, 
with one controlling substrate resistance and the other charac- 
terized by ambipolar diffusion. It was also found that steady- 
state funneling is more difficult to induce in the p + /n diode 
than in the n + /p diode. The carrier density exceeding the doping 
density in the substrate and at the DRB is not a sufficient 
condition to collapse a DR. A simple necessary condition for a DR 
collapse (or funneling) is derived in terms of ambipolar diffu- 
sion currents and is a statement regarding the spatial distribu- 
tion of carrier generation. The condition is satisfied if carrier 
generation is sufficiently close to the DR, but does not require 
generation inside of the DR. Quantitative predictions agree well 
with computer simulation results. 
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PREFACE (continued) 


This is the first rigorous (albeit steady-state) analysis of 
funneling in three dimensions, and may help to dispel some myths. 
Every point in the device lies on some equipotential surface, but 
a common misconception is that one such surface, called a fun- 
nel", is distinguished from the others by containing the region 
where there is a substrate electric field. In this picture, the 
electric field is in a region that extends a "funnel length" from 
the DRB into the substrate. In reality, the electric field is not 
confined to such a region and there is no unambiguous funnel. The 
region containing the strongest substrate electric field is 
typically adjacent to the electrode, where the carrier-density- 
modulated conductivity is smallest. This is seen under transient 
as well as steady-state conditions. The total substrate voltage 
drop measures the extent of DR collapse, but the distribution of 
this potential within the substrate merely responds to the 
carrier-density-modulated conductivity. Selection of a surface to 
be called a funnel is arbitrary, and the concept of a funnel was 
not found to be useful. Another common misconception is that 
funneling requires that carriers be generated inside the DR. In 
reality, carriers generated outside but close to the DR can also 
induce funneling. This is also seen under transient as well as 
steady-state conditions. 


The level of rigor accounts for the length of this analysis. 
Readers that are not interested in mathematical theory should be 
able to understand Chapters 1 and 6 without reading the other 
chapters . 
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1. INTRODUCTION 


This publication analyzes simple silicon diodes exposed to 
steady-state photon irradiation. Funneling (a term used by the 
single-event-effects community [1]) occurs when carriers are 
generated in sufficient quantity near a p-n junction depletion 
region (DR) that the DR becomes flooded and partially, or com- 
pletely, collapses. Some or nearly all voltage (including the 
built-in potential) normally across the DR is now across a sub- 
strate or epi layer, resulting in an electric field that enhances 
charge collection. This can occur under steady-state as well as 
transient conditions. The two types of conditions have some 
common qualitative characteristics, and concepts derived for the 
simpler steady-state problem can add physical insight into the 
more difficult transient problem. Theoretical transient models 
that exist at this time are unconvincing, and the primary motiva- 
tion for the present steady-state analysis is to obtain physical 
and mathematical guidance for a future transient analysis. There- 
fore special attention is given to the extremely high irradiation 
intensities needed to produce steady-state funneling, such as 
might be produced by a laser having a pulse width longer than the 
device relaxation time. The analysis is not limited to such high- 
intensity conditions, but these are the only conditions under 
which the conclusions derived here differ significantly from 
those derived from the classical theory. Even when classical 
theory is known to apply, the treatment of three-dimensional 
geometries presented here may be found to be useful. 

As shown in Figure 1.1, the simple silicon diode considered 
consists of a uniformly doped substrate between a p-n metallurgi- 
cal junction (MJ) and an ohmic contact (electrode) . The DR bound- 
ary (DRB) separates a strong space-charge region (the DR) from a 
quasi-neutral region. The simpler term "substrate" will refer to 
the quasi-neutral region from now on. Steady-state photogenera- 
tion occurs in the DR and/or substrate, and the generation rate 
density is assumed to be a known function (called the generation 
rate function) of the spatial coordinates. The figure shows an 
n + /p device, but results are also given for the p + /n device. The 
high-resistance region (HRR) , ambipolar region (AR) , and boundary 
(ARB) shown in the figure are discussed later. 

The nonlinear drift-diffusion equations are simplified by 
assuming constant mobilities in the substrate (although electric 
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Figure 1.1: Qualitative sketch of an n+/p diode showing a metal- 
lurgical junction (MJ) , a depletion region (DR) and its boundary 
(DRB) , an ambipolar region (AR) and its boundary (ARB) , and a 
high-resistance region (HRR) . The current I is positive when 
directed downward. 
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field dependent mobilities are used in the DR) and neglecting 
recombination (except at the electrode) . From this point on, the 
analysis is fairly rigorous. Compared to the classical analysis, 
the analysis given here is more general in one sense but more 
limited in another. It is limited to cases where recombination 
can be neglected. It is more general in the sense that it applies 
to a wide range of operating conditions including (but not limit- 
ed to) those that produce currents large enough for the classical 
law of the junction to break down, and that produce strong elec- 
tric fields in the substrate. Furthermore, the analysis applies 
to arbitrary substrate geometries and does not require that the 
DRB be a single connected surface. It can be the union of any 
number of disconnected surfaces (i.e., an array of DRBs) provid- 
ing that the same carrier density and potential boundary values 
are common to all surfaces. Similarly, the electrode can be the 
union of any number of disconnected ohmic contacts. However, if 
the DRB and/or electrode consist of several disconnected sec- 
tions, currents through the individual sections are not solved. 
Sums of currents (summed over the individual sections) are 
solved. 

The complete analysis consists of several distinct parts. One 
part, called the "DR analysis" solves the boundary value problem 
describing the DR. Another part, called the "substrate analysis" 
does the same thing for the substrate. The last part merely 
combines and solves the simultaneous equations provided by the 
other parts. Taken individually, the DR and substrate analysis 
are fairly general and can probably find applications in subjects 
other than an irradiated diode. 

The DR analysis was originally presented in a publication that 
few people know about [2]. The results as originally presented 
were so complex that they were virtually unusable. These results 
are greatly simplified in Appendix A, and apply to a broad range 
of conditions (high or low injection levels, with or without 
velocity saturation) . The substrate analysis (Chapters 3 and 4) 
applies to any substrate geometry and is more rigorous than 
analysis used in the past. It is never assumed in advance that 
one or another current component (electron or hole, drift or 
diffusion) can be neglected. It is sometimes concluded that one 
or another current component can be neglected, but the conclusion 
is derived (rather than assumed) and the conditions under which 
the conclusion is valid are quantified. Two special functions 
were found to be vital to the substrate analysis. These functions 


3 



are discussed extensively in Appendices B and C, which also 
contain subroutines for numerical evaluation. 

Solutions are expressed in terms of equilibrium resistance (the 
resistance between electrode and DRB that would occur if there 
were no excess carriers) , diffusion currents (predicted by the 
linear diffusion equation with simple boundary conditions) , and a 
nameless quantity derived from the photogeneration rate function. 
These quantities implicitly contain the required geometric data 
and substitute for physical dimensions in the formal solutions 
(e.g., instead of specifying a length and area, we specify an 
equilibrium resistance) . The advantage of this approach is that 
the equations are geometrically covariant, in the sense that the 
same equations are used for all geometries. Final numerical 
calculations are geometry specific and straightforward in one 
dimension. The three-dimensional case is made tractable by con- 
fining our attention to a special family of photogeneration rate 
functions, constructed so that all relevant functions of the 
spatial coordinates can be expressed as functions of a suitably 
chosen generalized coordinate (fitting is necessary if a given 
generation rate function does not belong to the family) . Some 
manipulations then show how numerical estimates can be obtained 
from the same calculations that would be used in one dimension. 
The user must provide an equilibrium resistance estimate and a 
fitting function representing photogeneration. All other calcula- 
tions, including diffusion current estimates, are first formally 
derived and then summarized in a "cookbook" recipe. 

The equations used for the substrate are familiar to everyone, 
and an earlier publication [2] provides the complete list of DR 
equations. After listing these equations, the analysis is mathe- 
matical. This explains the scarcity of references. Although this 
is a mathematical analysis and very few physical arguments are 
used in the derivations, physical interpretations are given for 
some of the mathematical results. No apology is given for the 
fact that the analysis is lengthy. This is unavoidable because we 
are solving a set of simultaneous nonlinear partial differential 
equations in three dimensions. The final result will be a numeri- 
cal algorithm for constructing the diode I-V curve corresponding 
to a given generation rate function. This can be done by a com- 
puter simulation, which can also treat diodes that are not sim- 
ple, and does not require the user to provide the resistance 
estimate and fitting function discussed earlier. The value of 
this analysis is physical insight, including verification of the 
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statements made in the paragraphs below. Numerical examples in 
the last chapter provide a visual illustration of predicted 
physical results. 

The present work finds that, when funneling is sufficiently 
strong, the ambipolar diffusion equation fails to provide a good 
approximation for the carrier density function, even when the 
predicted (via the ambipolar equation) carrier density is orders 
of magnitude greater than the doping density. The failure of this 
approximation is due to strong substrate electric fields. A more 
accurate equation is provided for quantitative calculations, but 
a simpler “generalized ambipolar approximation" is useful for 
visualization, and is described in the following way. 

The substrate divides into two subregions (see Figure 1.1). 
Adjacent to the electrode is an HRR characterized by a small 
excess carrier density and strong electric field. This region 
forms because funneling-induced substrate fields drive minority 
carriers up from the electrode. There are virtually no replace- 
ment carriers supplied by the electrode, so the region is deplet- 
ed of minority carriers. Quasi-neutrality insures that the region 
is also depleted of excess majority carriers. The conductivity is 
much less than in the high-density region above the HRR, so 
nearly all the substrate voltage drop is across the HRR. The 
region above the HRR is the AR and is characterized by a high 
carrier density and weak electric field. The ambipolar diffusion 
equation applies (approximately) to this region, but boundary 
conditions must be modified to account for the ARB that separates 
the AR from the HRR. It might be noted that the formation of an 
HRR and AR is very simple to derive in one dimension, if there is 
no photogeneration, and we assume in advance that the minority 
carrier current is negligible [3). The present work derives this 
result in three dimensions, with photogeneration, and without the 
up-front assumption. 

The HRR controls substrate resistance, while the ARB affects 
carrier density in the AR as if the electrode had been moved 
closer to the DRB. Furthermore, when funneling is sufficiently 
strong, the strong HRR electric field can drive nearly all minor- 
ity carriers to the DRB. Replacing the electrode with a high-low 
junction, which blocks the minority carrier current, will have 
little effect because this current is blocked anyway. The device 
is in saturation during sufficiently strong steady-state funnel- 
ing, i.e., nearly all liberated charge is collected. (This is one 
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distinction between the steady-state and transient cases. For the 
latter case, tunneling is strong during part of the charge col- 
lection time at most, and the collected charge can be less than 
the total amount liberated.) It should be noted that even when 
funneling is not strong enough to produce saturation, it can 
still be important enough make the device I-V curve significantly 
different than the classical prediction. 

The strong electric field in the HRR can affect mobility 
(velocity saturation) and it is reasonable to question the valid- 
ity of ignoring this effect and assuming constant mobilities in 
the substrate. It turns out that the currents are insensitive to 
this effect, because the field in the HRR does not become this 
strong until the device is well into saturation. The carrier 
density in the HRR does respond to this effect and the minority 
carrier density changes from one negligible value to some other 
negligible value. The good quantitative agreement between predic- 
tions given here and those given by a computer simulation that 
includes electric field dependent mobilities, indicates that it 
is not necessary to use electric field dependent mobilities in 
the substrate. 

A comparison is made between n + /p and p / n diodes having the 
same geometry and doping (except that n-type and p-type are 
interchanged), subject to the same bias voltage (except for a 
change in polarity) , and exposed to the same generation rate 
function. It was found that funneling is more difficult to induce 
in the n-type substrate. This observation goes beyond the simple 
fact that less mobile minority carriers are less responsive to a 
substrate electric field. In fact, the currents need not be 
greatly different and, depending on the bias voltage, either 
diode can have the larger current. The observation is that it+is 
more difficult to create a substrate electric field in the p /n 
device. In one numerical example, the voltage across the p-type 
substrate was 1.63 volts, compared to only 0.11 volt across the 
n— type substrate. The DR was greatly collapsed in the former 
case, but nearly intact in the latter case, even though the 
carrier density greatly exceeded the doping density in the sub- 
strate and at the DRB (implying that this is not a sufficient 
condition for a DR collapse) . A wide HRR occurred in the former 
case but not in the latter case; this compensated for the sub- 
strate voltage drops, so that the currents differed by less than 
22%. A simple necessary condition for saturation (or a DR col- 
lapse) is derived in terms of ambipolar diffusion currents, and 
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is consistent with the conclusion that tunneling is more diffi- 
cult to induce in the p + /n device. 

Readers that are not interested in mathematical theory can go 
directly to Chapter 6 beginning on page 73. 
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2. PRELIMINARY DISCUSSION AND GOVERNING EQUATIONS 


The analysis consists of several distinct steps. One step 
solves the equations describing the quasi-neutral region, which 
we have been calling the substrate. This is the region between 
the electrode (denoted S 1 for brevity) and the DRB (denoted S 2 ) • 
In this context, "solve" means that the electron and hole cur- 
rents are expressed in terms of the (unknown) carrier density and 
potential boundary values at S 2 . If we were treating a simple 
resistor, the equation V=IR with R known would be called the 
solution. The solution for the semiconductor substrate is worked 
out in Chapters 3 and 4 . Another step solves the equations de- 
scribing the DR. Again "solve" means that currents are expressed 
in terms of boundary values or vice-versa. This step was already 
done in a previous publication. The results were very messy and 
are simplified in Appendix A. The third and last step combines 
and solves the simultaneous equations for the currents and bound- 
ary values. This step is analogous to using Kirchhoff's laws to 
solve the problem of two resistors in series, and is worked out 
in Chapter 5. 

Because the DR analysis was already done, only the equations 
describing the quasi-neutral region need to be listed here. We 
start with the well-known equations which, under steady-state 
conditions with negligible recombination, reduce to 


J h ■ 9 D h [“ 9 rad P - (P + P Q ) grad U/V T ] (2.1a) 

J e = D e t9 rad N - (N + n Q ) grad U/V T ] (2.1b) 

div J h = q g (2.2a) 

div J e = - q g (2.2b) 

- e div grad U = q (P - N) (2.3) 

D e = V T ^e ' D h - V T **h (2.4) 
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where 


n , p = equilibrium electron and hole densities, respectively 
N, p = excess electron and hole densities, respectively 
D = diffusion constants for electrons and holes, respectively 

n h = mobilities for electrons and holes, respectively 
V^= thermal voltage (about 0.026 volts at room temperature) 
q = elementary charge 

j t = electron and hole current densities, respectively 
U = electric potential 
e = dielectric constant 
g = generation rate function 


The standard quasi-neutral approximation is obtained by regard- 
ing e as sufficiently small compared to other relevant constants 
that the solutions to the equations can be approximated by the 
solutions obtained in the limiting case as e approaches zero. In 
this limit, (2.3) is replaced with P=N and (2.1) and (2.2) are 
used to solve for* both P end U. 

Boundary conditions should also be stated. The reference poten- 
tial is chosen so that U=0 on S lt the semiconductor side of the 
electrode-semiconductor interface (contact potentials between 
electrodes and semiconductor will be included in Chapter 5) , 
where we also have P=0 . The values of P and U on S 2 are denoted 
P 2 and V 2 respectively, which are regarded as constants (in the 
spatial coordinates) on S 2 and represent some kind of spatial 
average on S 2 . All other boundary surfaces are insulated and 
assumed to be reflective for both electron and hole currents. 
This implies that the insulated boundaries are reflective for 
both P and U. 

Although not essential, it is notationally convenient to be 
definite as to whether the substrate is an n- or p-type. Only 
one case need be considered in detail because analogous results 
apply to the other case. All discussions and analysis will refer 
to the p-type substrate. Final equations will be listed for the 
n— type case in Sections 3.8 and 5.5. 

It is convenient to omit the equilibrium minority carrier 
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density n Q in (2.1b). This term produces a theoretically pre- 
dicted reverse current when the p-n junction is reverse-biased so 
that P 2 «-n Q . But this small current is not important because it 
is dominated by other currents (such as those associated with 
thermal generation/recombination in the DR) that are not included 
in this analysis. Therefore there is no compelling reason to keep 
the n Q and we will leave it out. 

The boundary value problem governing the p-type quasi-neutral 
region is now written as 


J h = 3 D h C“ ^rad P - (P + p 0 ) grad U/V T ] (2.5a) 

J e = q ° e [grad P - P grad U/V T ] (2.5b) 

div grad P + grad P • grad U/V T + (P + p Q ) div grad U/V T 


= “ g/D h (2.6a) 

div grad P - grad P • grad U/V T - P div grad U/V T 

= - g/D e (2.6b) 

P = 0, U = 0 on S ± (2.7a) 

P = P 2 , U = V 2 on S 2 (2.7b) 


grad P • n = 0, grad U • n = 0 on insulated boundaries (2.8) 


where n is the normal unit vector. The boundary value problem 
(2.5) through (2.8) is the mathematical definition of a "simple 
substrate" (for the p-type case under steady-state conditions) . 
Although a simple substrate can only approximate a real physical 
system (at best) , the equations themselves can be exactly solved 
for some special cases. An equation will be called exact if it is 
an exact mathematical result of these equations, regardless of 
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how well it represents a real physical system. 

The objective of the next two chapters is to solve these equa— 
tions so that the surface integrated currents are expressed in 
terms of P 2 and V 2 . Chapter 5 will solve for P 2 , V 2 , and all cur- 
rents. The surface integrated currents are defined by 


* h ,i s } s 'h • as - 

- q D h [grad P + (P + p Q ) grad U/V T ] • ds (i=l,2) (2.9a) 

J s i 


..i " [ s J e ' 

J s i 


ds = 


q D < 


[grad P — P grad U/Vip] • ds (i— 1,2) (2.9b) 


[ T,i _ I h,i + X e, i ^ lf2 ^ 


( 2 . 9c) 


where the unit normal vector in all surface integrals is an outer 
normal, i.e. , directed away from the substrate interior. A sur 
face integrated current is positive if positive charge moves 
toward the surface from the substrate interior. 

Adding (2.5a) to (2.5b) and adding (2.6a) to (2.6b) produces a 
result that can be written as 


J T = - a grad U H 


( 2 . 10 ) 


div J T = 0 


( 2 . 11 ) 


where 


J T s J e + J h 


( 2 . 12 ) 
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is the total current density and 


a = q [M e P + M h (P + P G )] = (q/V T ) (D g + D h ) (P + A Q ) (2.13) 

is the conductivity, with the constant A Q defined by 

A o = D h Po/( D e + D h> * ( 2 . i4 ) 

The function U H is defined by 

U H ■ U - (2V T /p 0 ) (p 0 /2 - A q ) ln(l + P/A q ) . (2.15) 


Note that (2.10) and (2.11) are simply Ohm's law except that the 
"potential" is U H instead of the actual potential U. The inte- 
grated form of Ohm's law is V=IR or 


V 2 - (2V t / Pq ) ( p 0 /2 - A 0 ) In (1 + P 2 /A q ) = I Tfl R (2.16) 

where R is the resistance between S-^ and S 2 produced by the 
conductivity a. This equation has limited computational applica- 
tions because the carrier-density-modulated resistance R is un- 
known. The equation does have some applications, which will help 
to reach some conclusions in Sections 3.2 and 3.4. 

Some constants and functions are defined below for later use. 

The equilibrium conductivity ct_. and ambipolar diffusion coeffi- 
. * , U 
cient D are defined by 


a o - <3 ^h Po = «3 / v T> D h Po 


(2.17) 


1/D* = (1/D h + l/D e )/2 . 


(2.18) 
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The unit function fi u and the function 0 are defined by the 
ary value problems 

div grad n u = 0 in substrate 
n u = 0 on s 1 

n u = 1 on s 2 

grad n u • n = 0 on insulated boundaries 


div grad 0 = - g/D* in substrate 
<p = 0 on Sj 

0=0 on S 2 

grad 0 • n = 0 on insulated boundaries . 

Associated with these functions are the parameters R Q , G ±l 
defined by 


1/R 0 s - 


grad n u • ds = ct ( 


grad n u • as 


= - q D 


grad 0 • ds (i=l,2) 


The G's are related by 


bound- 

(2.19a) 

(2.19b) 

(2.19c) 

( 2 . 19d) 

(2.20a) 
(2.20b) 
(2.20C) 
(2 . 20d) 
and G 2 

( 2 . 21 ) 

(2.22a) 
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Gi + G 2 = q 


g d 3 x . 


(2.22b) 


J sub 


Each of these parameters has a physical interpretation. R Q is the 
electrical resistance between S.^ and S 2 produced by the uniform 
equilibrium conductivity a Q . G^ (i=l,2) is the absolute value of 
the ambipolar diffusion current through that would occur if 
the carrier density satisfied the ambipolar diffusion equation 
with S 1 and S 2 both acting as sinks for excess carriers (i.e., if 
P =<p) . These parameters are constants in the sense that they do 
not depend on spatial coordinates or on the boundary values P 2 or 
V 2 . However, they do depend on operating conditions. In addition 
to the obvious dependence that G^ has on g, there is also an 
implicit dependence due to the fact that the location of the 
boundary S 2 , which defines the geometry, can vary due to varia- 
tions in the DR width. It will not be necessary to consider 
variations in the boundary S 2 until we get to Chapter 5. Chapters 
3 and 4 will proceed as if the boundary location and boundary 
values are known and fixed. The parameters R Q and the G's are 
regarded as known when the boundary location is given. Chapter 4 
will show how the G's can be calculated from a particular type of 
function used to fit g. 
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3. SUBSTRATE ANALYSIS: A SPECIAL CASE 


3 . 1 Introduction 

We begin with a practice problem in which there is no photogen— 
e *’ a ^i° n in the substrate. Although simpler than the more general 
case, this special case is far from trivial because carriers can 
be injected through S 2 * S 2 will be a p- n junction DRB in Chapter 
5, but can presently be the boundary of any physical structure, 
because the boundary values P 2 and V 2 are arbitrary. In particu- 
l ar f it can represent a high-low junction, a forward biased p-n 
junction injecting minority carriers into the substrate, or a 
reverse-biased p-n junction injecting majority carriers into the 
substrate via photogeneration within the DR. Some concepts ap- 
plicable to more general conditions are most easily discovered by 
starting with this problem, because the analysis is not burdened 
by a lot of mathematical complexity and an exact solution can be 
found. Of special interest is the formation of an HRR and AR 
(discussed later) when V 2 is large and positive (a p-type sub- 
strate is assumed) . This situation (funneling) occurs if carriers 
are generated within a reverse-biased DR fast enough to flood it, 
causing it to collapse so that much of the applied plus built-in 
voltage is across the substrate. 

The analysis to follow regards the location of S 2 and the 
boundary values P 2 and V 2 as given constants. The equilibrium 
resistance R Q is regarded as known, so the currents are consid- 
ered to be solved when expressed in terms of P 2 , V 2 , and R Q . 


3 . 2 Solution for P and U 


By adding (2.6a) to (2.6b) while using g=0, we obtain 


div grad [P + (p Q /2V T ) U] = 0 . 


(3.1) 


Comparing the boundary value problem satisfied by the expression 
in brackets to (2.19), we find that 


P + (p Q / 2 V T ) U = n 
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where 

n s [p 2 + (p g / 2V t ) v 2 ] n u 


is regarded as a known function of the spatial coordinates. Using 
(3.2) to eliminate U in either (2.6a) or (2.6b) gives 


div [(P + P c /2) grad (P - n) ] = 0 


(3.4) 


The solution to this equation is P satisfying 


P + ( Po /2 - A) ln(l + P/A) = n 


where A is a constant. Substituting (3.5) into (3.4) verifies 
that (35) is a solution. The boundary conditions are satisfie 
at Sl . The constant A is selected so that the boundary conditions 
are also satisfied at S 2 . Evaluating (3.5) at S 2 , we fin a 
satisfies 


(p 0 /2 - A) In ( 1 + P 2 /A) = (P q / 2V t ) V 2 


(3.6) 


and can be calculated from either 


A = (p 0 /2) [1 - (V 2 /V T ) E] if V 2 * - 2 V t P 2 /p 0 (3-7a) 


or 


p (gl/E — 1) i if V 2 =$= — 2Vip ^ 2^0 


(3.7b) 


where 


E = [ln(l + 2P 2 /p 0 )] _1 if V 2 = 0 and P 2 > 0 


(3.8a) 
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E ■ H(Z 1 ,Z 2 ) if V 2 0 and V 2 =}= - 2V T P 2 /P 0 (3.8b) 


Z'l = (V t /V 2 ) (1 + 2P 2 /p 0 ) , Z 2 = V T /V 2 


(3.8c) 


and the special function H is defined by 


H(Z 1 ,Z 2 ) = E if and only if exp(l/E) = (E-Z-jJ / (E-Z 2 ) . (3.9) 

Equations (3.7a) and (3.7b) give the same result in theory, but 
(3.7b) should be used if (V 2 /V T )E is so nearly equal to 1 that 
(3.7a) requires more numerical precision than is -available. 
Otherwise, (3.7a) can be used. 

Properties of the function H are discussed in Appendix B, which 
also contains a subroutine for numerical evaluation. Although not 
obvious from a casual inspection of (3.9), there is a problem if 
1+Z 1 -Z 2 =0. As 1+Z 1 ~Z 2 approaches zero, H(Z lf Z 2 ) becomes posi- 
tively or negatively infinite, depending on whether the approach 
is from above or below. This problem case occurs when 
p 2 + (Po/ 2V t) v 2 =0 so that n=0. The solution given by (3.5) does not 
apply to this case and must be replaced with 


(P + p 0 /2 ) 2 = [(P 2 + p 0 /2 ) 2 - ( p 0 / 2 ) 2 ) n u + (p 0 / 2 ) 2 if n=o 

which is easily verified by substituting it into (3.4). This 
problem case will occur if S 2 is an electrode (P 2 =0) and shorted 
to (V 2 =0) . But even if there is photogeneration in the sub- 
strate, this case is still not very interesting because, accord- 
ing to (2.16), the terminal current is zero. Other than this 
uninteresting example, the problem case would be associated (at 
least in concept) with a forward-biased DR (V 2 <0) with the for- 
ward biasing strong enough to produce a large voltage drop across 
the highly conductive substrate. The current would quickly de- 
stroy the device. The problem case is not expected in applica- 
tions of interest, so we will always use the solution given by 
(3.5) with A solved from (3.7). 
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P is solved from (3.5) and U is solved from (3.2). The solution 
for P can be written more explicitly by defining another special 
function F by 


F(X 1 ,X 2 )= Y if and only if Y + (1-X 1 ) ln(l + Y/X 1 ) - X 2 . (3.10) 


Properties of F are discussed in Appendix C, which also contains 
a subroutine for numerical evaluation. Comparing (3.5) and 
(3.10) , we get 


P = (p 0 /2) F(2A/p 0 , 2n/p 0 ) . (3.11) 


3.3 Solution for the Currents 

By taking the gradients of (3.2) and (3.5) and combining equa- 
tions we get 

grad P = [(P+A)/(P + P 0 /2)] grad SI (3.12) 

grad U = (2V T /p Q ) [ (p c /2 - A)/(P + p 0 /2)] grad n . (3.13) 

Substituting these gradients into (2.9) gives 


I h ,2 = - J h,l = 2 <J D h I 1 - A/p o> 
I«,2 - - r e,l “ - 2 « D e < A /Po> 


grad Si • ds 


grad H * ds 


JS 


1 


and combining with (3.3) and (2.21) gives 
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J h , 2 = - J h,l = - (1 " A/ Po ) (V 2 + 2 V t P 2 /p 0 )/R 0 (3.14a) 

J e,2 = " I e, 1 = (°e/ D h) ( A /Po> ( V 2 + 2V T P 2 /P 0 ) / R o (3- 14b) 

which, together with (3.7), completes the solution for the cur- 
rents . 


3 . 4 The Nominal Ambipolar Approximation 

The behavior of P is easier to visualize if transcendental 
equation (3.5) is approximated by a simpler equation. The sim- 
plest approximation, which has some applications when P 2 >>p Q , is 
the nominal ambipolar approximation obtained by neglecting U in 
(3.1) to get 


P « P 


(3.15) 


where P* is defined (when g=0) by the boundary value problem 


div grad P = 0 in substrate 


(3.16a) 


P = 0 on S, 


(3.16b) 


P = P 2 on S 2 


(3.16c) 


grad P • n = 0 on insulated boundaries (3.16d) 


Comparing (3.16) and (2.19), we find that 


21 



(3.17) 


= p 0 n 


u 


We can use (3.2) and (3.3) to conclude that the nominal ambipolar 
approximation (3.15) is valid if the ambipolar condition 


P 2 » (p 0 /2V T ) | V 2 | (ambipolar condition) (3.18) 


is satisfied. 

Some of the older literature gives a misleading impression 
regarding ambipolar diffusion. The impression given is that 
electrons and holes interact so strongly, through their mutual 
attraction, that they move together and do not respond to applied 
fields. This picture accounts for U being absent in the equation 
governing P, but also predicts that J T =0 (because electrons and 
holes move together) . The assertion Jip=0 has also been supported 
by analysis of a strongly symmetric problem (cylindrical symmetry 
with no longitudinal flow) . But such strong symmetry has some 
properties (e.g., the divergence of a bounded vector field 
uniquely determines the vector field) that do not apply to more 
general cases. The conclusion does not apply if the symmetry is 
weaker (e.g., cylindrical symmetry but with longitudinal flow) or 
if there is no symmetry. In the more general case, electrons and 
holes can move very differently from each other while maintaining 
quasi-neutrality, if carriers moving out of a volume element are 
replaced by others moving in. While it is true that the carrier 
density function is insensitive to weak applied fields, carrier 
motion is very responsive. This response can be seen from (2.16). 
R is insensitive to V 2 , so the total current is nearly linear in 
V 2 - Even when the ambipolar approximation is known to apply, we 
should avoid additional approximations derived from the idea that 
electrons and holes move together and independently of applied 
fields. 


3.5 A Generalized Ambipolar Approximation 

It is possible to modify the nominal ambipolar approximation to 
include some cases violating the ambipolar condition (3.18). We 
do assume throughout this discussion that P 2 >>p Q /2. There are 
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four cases that can be considered. For the first case, V 2 is 
positive but small, where "small" means several times V T . For the 
second case, V 2 is negative but small in absolute value. The 
nominal ambipolar approximation should apply to both of these 
cases. For the third case, V 2 is negative but large (>>V T ) in 
absolute value. This case is not of practical interest. A physi- 
cal arrangement producing this case is one in which S 2 represents 
a forward-biased p-n junction with an applied voltage strong 
enough to produce a large ohmic voltage drop across the highly 
conductive substrate. The large currents will quickly destroy the 
device. For the fourth case, V 2 is positive and large. This case 
can occur without destroying the device because a current limit- 
ing HRR forms (discussed below) . A number of physical arrange- 
ments can produce the fourth case. Of special interest here is 
the one in which S 2 represents a reverse-biased p-n junction with 
photogeneration within the DR strong enough to collapse it, so 
that much of the applied plus built-in voltage is across the 
substrate (funneling) . Given that P 2 >>p Q /2, the fourth case is 
the only case of practical interest where the nominal ambipolar 
approximation fails. The objective of this section is to general- 
ize the ambipolar approximation to include this case. The remain- 
der of this section assumes that V 2 is positive. 

An approximation for P can be derived by taking the gradient of 
(3.5) to get 


grad P = [ (P + A)/(P + p 0 /2) ] grad n . (3.19) 

It can be shown that a positive V 2 implies that A satisfying 
(3.6) also satisfies 


0 < A < p 0 /2 if V 2 > 0 


(3.20) 


By assumption, P 2 >>p 0 /2. Therefore there is some region adjacent 
to S 2 where P>>p Q /2 and P>>A, so that the bracket in (3.19) is 
nearly unity, i.e., gradPssgradfl , implying that P and n differ (in 
this region) by an additive constant. The additive constant can 
be evaluated by noting that the region includes S 2 . The result is 
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P » n - (P c / 2 V t ) v 2 = [P 2 + (Po/ 2V T> V 2 1 n u ■ (Pq/ 2V t) V 2 * 


This equation is valid in a region sufficiently close to S 2 to 
satisfy P>>p Q /2. Any points where the right side of the above 
equation is zero cannot be in this region. The ARB is mathemati- 
cally defined to be the set of points where the right side of the 
above equation is zero, i.e., the constant n u surface character- 
ized by 


n u = (P 0 /2V t ) V 2 /[P 2 + (P q / 2V t ) V 2 ] defines ARB . (3.21) 

The AR is mathematically defined to be the region between the ARB 
and S 2 . Excluding a transitional region adjacent to the ARB, the 
AR is characterized by P>>p 0 /2 so that 


P » [P 2 + (P 0 /2V t ) V 2 ] Il u - (P 0 /2V T ) V 2 in AR . (3.22a) 

The HRR is mathematically defined to be the region between the 
ARB and the electrode S^. It can be shown from the exact equa- 
tions that, excluding a transitional region adjacent to the ARB 
(where P can be several times p^/2) , the HRR is characterized by 
P«p 0 so that 


P a 0 in HRR . 


(3.22b) 


The HRR is characterized by a low conductivity («a Q , which is 
small compared to the conductivity in the AR) and a large (nearly 
all of V 2 ) potential drop when V 2 »V T . This motivated the name 
"high-resistance region". This region limits the current so that 
a large V 2 can occur without destroying the device. The AR region 
is characterized by a high conductivity and small (several times 
V T ) potential drop. These are the conditions appropriate for 
ambipolar diffusion and motivated the name "ambipolar region" . 

We temporarily drop the assumption that V 2 is positive and 
define the generalized ambipolar approximation to be (3.22) when 
V 2 is positive and (3.15) otherwise. Reinstating the assumption 
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that V 2 is positive (so that the ARB exists), it is evident from 
(3.21) that the ARB becomes in the limit of small V 2 . In this 
same limit, the generalized approximation (3.22) reduces to the 
nominal approximation (3.15). 

There is a physical explanation for the absence of excess 
carriers in the HRR. An HRR with sufficient width to be depleted 
of excess carriers (i.e., the HRR is distinguishable from the 
transitional region) forms when V 2 is large enough for the gener- 
alized ambipolar approximation to significantly differ from the 
nominal ambipolar approximation. But electric fields strong 
enough to make the nominal approximation fail are also strong 
enough to push electrons away from the electrode. The electrode 
supplies virtually no electrons, so there is a region near the 
electrode that is virtually depleted of electrons. Quasi-neutral- 
ity implies that this region is also virtually depleted of excess 
holes . 

An alternate definition for the ARB, mathematically equivalent 
to (3.21), can be stated in terms of the slope of P. This alter- 
nate definition makes the ARB easier to visualize. The general- 
ized and nominal ambipolar approximations predict the slope of P 
near S 2 to be given by 


grad P « [P 2 + (p Q /2V T ) V 2 ] grad n u (generalized) (3.23a) 
grad P « P 2 grad fl u (nominal) (3.23b) 


so that the generalized approximation predicts a steeper slope 
than the nominal approximation. The nominal approximation can be 
modified to give the generalized approximation by moving the sink 
boundary from the electrode to the ARB. Moving the sink boundary 
closer to S 2 produces a steeper slope. The ARB can be visualized 
(and defined) as the location where the sink must be placed to 
produce the correct (steeper) slope. 

The generalized ambipolar approximation must be used with 
caution and should not be used in calculations that subtract 
nearly equal quantities and require high accuracy. For example, U 
is solved from (3.2) after P has been solved, but the exact 
solution must be used. Using the approximation for P will pre- 
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diet a zero electric field in the AR. This is not a good estimate 
of the electric field. The electric field is small in the AR only 
because the conductivity is correspondingly large, so even a 
small electric field is important and cannot be neglected. The 
generalized ambipolar approximation is an approximation for (3.2) 
and might be used instead of (3.2) (requiring that U be solved 
some other way), but cannot be used with (3.2). The approximation 
is useful for visualization, for predicting distinct regions 
where P has different behaviors, and for defining the ARB which 
separates these regions. But (3.11) is recommended for numerical 

calculations . 

The final observation made here concerns the electron current. 
The conditions (large V 2 ) that result in the generalized ambipo- 
lar approximation being significantly different than the nominal 
approximation also result in A being extremely small. The elec- 
tron current given by (3.14b) is extremely small. The physical 
explanation is the same as that given for the absence of excess 
carriers in the HRR. An electric field strong enough to cause the 
nominal approximation to fail is also strong enough to prevent 
electrons from reaching the electrode, so I e ,l~°* Th * s phy s i cal 
explanation also applies to the g+o case considered in Chapter 4. 
It is interesting to note that under large V 2 conditions, it 
makes no difference whether S-j^ is an electrode or a high-low 
junction that blocks the electron current because I &1 is virtu- 
ally zero anyway. 


3.6 Low-Ini ection-Level Conditions 

Low-injection- level conditions (LILC) occur when P«P G through- 
out the substrate. It is commonly assumed that LILC implies that 
the minority carrier diffusion equation (MCDE) gives a good 
approximation for P. It is interesting to determine whether this 
assumption is valid. It turns out that the assumption is invalid, 
but can still be used for the purpose of estimating total cur- 
rent. The meaning of this statement is explained below. It is 
also shown that the MCDE applies if and only if A»P. 

Given LILC, a necessary condition for the MCDE to apply can be 
determined by comparing the MCDE-predicted gradients of P at S x 
and S 2 to the actual gradients. The solution to the MCDE, for 
steady-state conditions with negligible recombination/generation, 


26 



is P 2 fl u (the same as the nominal ambipolar approximation) . The 
predicted gradient of P at either boundary is P 2 gradn u . The 
actual gradient is given by (3.19). Using (3.3) gives 


grad P = (2/p Q ) A [P 2 + (p D /2V T ) V 2 ] grad n u at S 1 

grad P = [ (P 2 + A)/(P 2 + p 0 /2)J [P 2 + (p 0 /2V T ) V 2 ] grad n u 

« ( 2/p 0 ) (P 2 + A) [P 2 + (p 0 /2V T ) V 2 ] grad n u at S 2 . 


One necessary condition for both of the above gradients to ap- 
proximately equal P 2 gradn u is A»P 2 , so that the coefficients on 
the two right sides will be nearly equal to each other. Another 
necessary condition is 


(2/ Po ) A [P 2 + (p 0 /2V T ) V 2 ] * P 2 
or 


2 P 2 / Po + V 2 /V T « P 2 /A . 


But P 2 /p 0 «l and P 2 /A«l, so |V 2 /V t |« 1. We conclude that LILC 
are not sufficient for the MODE to apply. It is also required 
that |V 2 |«V T . 

A different line of reasoning will conclude that, given LILC, 
we can pretend that the MCDE applies, even if it really does not, 
providing that our interest is in total current. Given LILC, the 
minority carrier drift current is negligible compared to the 
majority carrier drift current. If the diffusion currents are on 
the order of, or larger than, the majority carrier drift current, 
then the minority carrier drift current is negligible compared to 
all other currents, and the MCDE applies (implying that |V 2 |<<V^,, 
which is consistent with the statement that majority carrier 
drift is as small as diffusion) . If the diffusion currents are 
much smaller than the majority carrier drift current, then the 
MCDE does not apply. But we can pretend that it does, because 
nearly all current is majority carrier drift and error in the 
calculated diffusion current does not matter. Note that the MCDE 
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implies that A>>P. Therefore, when calculating total current with 
LILC given, we can assume that A>>P, even though the assumption 
may be wrong. Conditions under which the assumption is wrong are 
also conditions under which error in the assumption does not 
matter . 

It was shown above that the MCDE implies that A>>P. It is 
interesting that the implication also goes in the other direc- 
tion. Given that A»P, we can expand the logarithms in (3.5) and 
(3.6) so that the equations reduce to 


P + (p 0 /2 - A) (P/A) « n = [P 2 + (p 0 /2V T ) v 2 ] n u 
( Po /2 - A) (P 2 /A) * (P g / 2V t ) V 2 


and combining equations to eliminate A gives P«P 2 fi u . 


3.7 Summary of Results for the p-Type Substrat e 

The results are now summarized for the p-type substrate. The 
equilibrium conductivity a Q is qM^Po or ( < 3/ v T^ D hPo where the 
equilibrium hole density p Q can be equated to the doping density. 
The equilibrium resistance R Q is the electrical resistance 
between S-,^ and S 2 produced by the equilibrium conductivity. The 
constant A is calculated from either 


A = (p 0 /2) [1 - (V 2 /V T ) E] if V 2 + - 2 V t P 2 /P ( 


or 


A = P 2 (e 


1/E _ 


1 ) 


-1 


if V, - 2 V t P 2 /p c 


where 


E ■ [ln(l + 2P 2 /p 0 )] -1 if V 2 = 0 and P 2 > 0 
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E = H(Z 1 ,Z 2 ) 


if V 2 * 0 and V 2 + - 2 V t P 2 /p q 


Zi s (V t /V 2 ) (1 + 2P 2 /p 0 ) , Z 2 s V T /V 2 

and the special function H is defined in Appendix B. The two 
equations for A give the same result in theory, but the second 
should be used if (V 2 /V T )E is so nearly equal to 1 that the first 
requires more numerical precision than is available. Otherwise, 
the first can be used. The exceptional case where A and E are 
undefined is mathematically possible but should not be encoun- 
tered in practical applications. The currents are calculated from 


J h, 2 = - T h,l = " (1 " A/ Po ) (V 2 + 2 V t P 2 /P 0 )/R q 

X e,2 " " ^ e, 1 = (°e/ D h) ( A /P 0 ) (V 2 + 2V T P 2 /p 0 ) /R Q . 

The above equations complete the solution for the substrate in 
the case where there is no photogeneration in the quasi— neutral 
region. But it is interesting to also look at the function P. The 
exact solution is given by either 


P + (p 0 /2 - A) ln(l + P/A) = n 
or 

P = ( Po /2) F ( 2 A/p 0 , 2n/p 0 ) 

where 

n - [ p 2 + (Pq/ 2V t) V 2 1 n u 

with the unit function fl u defined by (2.19) and the special 
function F discussed in Appendix C. Approximations are available 
for P. First assume that P 2 <<p Q /2. Then either majority carrier 
drift is the dominant current, or P«P 2 n u and A>>P. Now assume 
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that P 2 »p 0 /2. For cases of practical interest such that V 2 <0, 
the approximation is 


P * p 2 n u if V 2 < 0 and P 2 » P 0 / 2 * 

If v 2 >0, an approximation is obtained by defining the ARB to be 
the constant fl surface characterized by 


n = (P d / 2V t ) v 2 . 


The AR is the region between the ARB and S 2 while the HRR is the 
region between the ARB and The approximation is 


P « n - (p 0 /2V T ) V 2 in AR if V 2 >0 and P 2 »p Q /2 


P « 0 in HRR if V 2 >0 and P 2 » P 0 /2 . 


The approximation is useful for visualization, but the solution 
in terms of F is recommended for numerical calculations. 


3.8 Analogous Results for t he n-Type Substrat e 

The analogous results are summarized for the n-type substrate. 
The equilibrium conductivity o Q is qji e n 0 or (q/V T )D e n Q where the 
equilibrium electron density n Q can be equated to the doping 
density. The equilibrium resistance R Q is the electrical resis- 
tance between and S 2 produced by the equilibrium conductivity. 
The constant A is calculated from either 


A = (n 0 /2) [1 + (V 2 /V T ) E] if V 2 + 2V T P 2 /n Q 


or 
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A = P 2 (e 1 / 2 - l)" 1 if V 2 =|= 2 V t P 2 /n 0 

where 

E ■ [ In (1 + 2P 2 /n 0 )]“ 1 if V 2 = 0 and P 2 > 0 
E ■ H(Z 1# Z 2 ) if v 2 ={= 0 and V 2 2V T P 2 /n Q 
Zj. = - (V T /V 2 ) (1 + 2P 2 /n Q ) , Z 2 - - V T /V 2 


and the special function H is defined in Appendix B. The two 
equations for A give the same result in theory, but the second 
should be used if (V 2 /V T )E is so nearly equal to -1 that the 
first requires more numerical precision than is available. Other- 
wise, the first can be used. The exceptional case where A and E 
are undefined is mathematically possible but should not be en- 
countered in practical applications. The currents are calculated 
from 


J h,2 - - J h,l = (°h/ D e) ( A /n 0 )(V 2 - 2 V t P 2 /n 0 )/R 0 

I e,2 = ” I e,l = ” (1 ' A / n o> < V 2 " 2V T p 2 / n o) / R o 
The exact solution for P is given by either 

P + (n 0 /2 - A) In ( 1 + P/A) = n 
or 

P = (n 0 /2) F(2A/n 0 , 2n/n 0 ) 

where 

n = [P 2 - (n D /2v T ) V 2 ) n u 
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with the unit function n u defined by (2.19) and the special 
function F discussed in Appendix C. Approximations are available 
for P. First assume that P 2 «n D /2. Then either majority carrier 
drift is the dominant current, or P«P 2 fi u and A>>P. Now assume 
that P 2 >>n 0 /2. For cases of practical interest such that V 2 >0, 
the approximation is 


p ~ p 2 n u if V 2 > 0 and P 2 » n Q /2 . 


If v 2 <0, an approximation is obtained by defining the ARB to be 
the constant n surface characterized by 


n = - (n D /2V T ) V 2 . 


The AR is the region between the ARB and S 2 while the HRR is the 
region between the ARB and S x . The approximation is 


P * fl + (n Q /2V T ) V 2 in AR if V 2 <0 and P 2 »n Q /2 

P » 0 in HRR if V 2 <0 and P 2 » n Q /2 . 

The approximation is useful for visualization, but the solution 
in terms of F is recommended for numerical calculations. 
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4. SUBSTRATE ANALYSIS: THE GENERAL CASE 


4 . 1 Introduction 


We now consider the general case in which there is photogenera- 
tion in the substrate. Unlike the special case in Chapter 3, 
exact solutions are not available for the general case. An exact 
analysis is used in Section 4.2 to express all currents in terms 
of I e 1 (a p-type substrate is assumed here) . Another exact 
analysis in Section 4 . 3 expresses I e 1 in terms of a new unknown 
function r, which will eventually be approximated. Function r is 
constructed in such a way that an estimate of I e 1 is insensitive 
to error in r. Section 4.4 gives an approximation for P, which is 
first used to approximate r, then I e 1 , and then the other cur- 
rents. A mathematical theorem in Section 4.5, a suitable restric- 
tion on g discussed in Section 4.6, and a numerical integration 
discussed in Section 4.7 make the approximations computationally 
manageable. Unlike Chapter 3, this chapter does not end with 
summary sections, because the final equations (including those 
for the n-type substrate) are summarized in Sections 5.3 and 5.5. 


4 . 2 Expressing Currents in Terms of Ie,l 

By adding (2.6a) and (2.6b) while using (2.18), we obtain 

div grad [P + (p Q /2V T ) U] = - g/D* . (4.1) 

Comparing the boundary value problem satisfied by the expression 
in brackets to (2.19) and (2.20), we find that 


P + (P q / 2V t ) U = n + <f> (4.2) 

where 

n s [P 2 + (p q / 2V t) V 2 J n u * (4.3) 


The two divergence equations (2.2a) and (2.2b) allow S 2 currents 
to be related to currents according to 
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c h,2 - 9 


sub 


g d 3 x - I h ,l = G 1 + G 2 ' *!>, 1 


x e,2 - - <J 


sub 


g d 3 x - I 6fl = - - G 2 - I e#1 . (4.4a) 


Taking the gradient of (4.2) and using (2.9) allows the Sj^ cur- 
rents to be expressed in terms of gradU and grad(n+0), which then 
allows I h to be expressed in terms of I e ^ as 


J h, 1 “ ( D h/ D e^ I e,l 2 q D h 


grad (n + 0) • ds 


and using (4.3), (2.21), and (2.22) gives 


I h / i = (V 2 + 2 V t P 2 /P 0 )/ r o + (! + D h/°e) G 1 + ^ D h/ D e) I e,l 


and the equation for Ijj 2 becomes 


J h,2 “ g 2 - < D h/ D e> G 1 

— (v 2 + 2 V(j P 2 /Po)/ R o ~ (^h/^e^ ^e,l 


(4.4b) 


4.3 Expressing Ie.l in Terms of r 

Using (4.2) to eliminate U in (2.6b) and rearranging terms 
gives 

div {(P + A) grad [ ) } = (A Q - A) div grad <p (4.5) 
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where 


[ ] = P + (P 0 / 2 " A) In ( 1 + P/A) - (n + $) 


and A is given by (3.6) . Note that A could have been replaced by 
other constants in the above equations. The motivation for the 
particular choice A will be clear later. We now define a new 
unknown function T by the boundary value problem 


div [ (P + A) grad T] = 0 in substrate (4.6a) 

T = 0 on S 1 (4.6b) 

T = 1 on S 2 . (4.6c) 

The present objective is to express I e x in terms of r, so that 
an approximation for r, which will come later, produces an ap- 
proximation for I e ^ . The divergence theorem together with (4.5) 
and (4.6) gives 


(1 - T) (P + A) grad [ ] • ds + 


[ ] (P + A) grad r • ds 


+ (A - A q ) 


(1 - T) grad <p • ds = (A Q - A) 


grad r • grad <p d 3 x 


where the surface integrals are on both S x and S 2 , and the volume 
integral is over the substrate. Using 

(P + A) grad [ ] = (P + p 0 /2) grad P - (P + A) grad (n + <p) 

together with (4.6b) and (4.6c) gives 
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(Po/2) 


grad P • ds = A 


grad (n +$) • ds 


] (P+A) grad r • ds 


+ (A 0 



“ 



A) 


grad <p • ds + 

grad r • grad <p d 2 x 

• 

S 1 

sub 


which expresses the left side in terms of known quantities and 
the unknown r. The motivation for selecting A instead of some 
other constant is that [ ]=0 on S x and S 2 . The equation reduces 
to 


(P 0 /2) grad P • ds 

K 


grad 


(n + <p) 


ds 


+ ( A q - A) 



p 

grad <p • ds + 

grad r • grad <p d 3 x 

\s 1 J 

sub 


so the unknown r appears only in a weight factor in a weighted 
average. This observation will be used in the next section, which 
produces an approximation for I e ,l* The above eg uat i° n can be 
expressed in terms of I e ^ using (2.9) and (4.2) with the result 


[Po /< 2 D eH l e ,l = A 


grad (n + <fi) • ds 


+ (A 0 - A) 


grad <p • ds + 


sub 


grad r • grad <f> d 3 x 


(4.7) 


36 



4 . 4 An Approximation for P and the Currents 


The role that r plays in (4.7) is most visible when the equa- 
tion is written in one dimension as 


[P 0 /(2 q D e ) ] I e fl (per unit area) = - A d(n+0)/dx Q 


+ (A 0 - A) 


- d0/dx o + 


(dr/dx) (d0/dx) dx 


(1 dim.) (4.8) 


where is at x=0, S 2 is at x=L, and d/dx Q is abbreviated nota- 
tion for the derivative evaluated at x=0. The normalization 
condition (4.6b) and (4.6c) can be written as 


*L 

(dr/dx) dx = 1 

J 0 


so dr/dx in (4.8) is the weight factor in a weighted average of 
d0/dx. Integrating (4.6) gives an alternate expression for the 
weight factor 


dr/dx = 


*L 

1/(P + A) dx 

0 


-1 

[1/(P + A) ] 


(1 dim.) 


(4.9) 


If V 2 is positive and large, A is very small and the weight 
factor is concentrated near x=0, where P=0. The weighted average 
reduces to the endpoint value at x=0 and I e is small. This is 
the expected result when V 2 is large. 

Weighted averages are usually insensitive to small errors in 
the weight factor, and this suggests that I e 1 can be approximat- 
ed by replacing the unknown r in (4.7) with an approximation. An 
approximation for r is obtained by replacing the unknown P in 
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(4.6) with an approximation. For LILC, we can assume that P<<A, 
so it does not matter how we approximate P in (4.6), as long as 
the approximation is consistent with P<<A. We therefore look for 
an approximation applicable to high-injection-level conditions 
(HILC) . 

The present objective is to find an approximation for P sp - 
plicable when P»p Q /2 throughout most of the substrate. A tenta- 
tive approximation is P T which is defined by 


P T + (p 0 /2 - A) In ( 1 + P T /A) = n + <p (4.10) 


and satisfies the required boundary conditions. To establish the 
credibility of the approximation P T , note that (4.5) can be 
manipulated into 


div [ (P + p 0 /2 ) grad P] = div [ (P + A Q ) grad (n + <p) ] (4.11a) 


while (4.10) can be used to show that 


div [ (P T + P 0 /2) grad P T ] = div [ (P T + A) grad (n + <t>) ] . (4.11b) 

The two equations differ only in that one contains A Q while the 
other contains A. The constant A Q is on the order of p Q /2, while 
A will be of the same order or smaller. For HILC, we will have 
P>>A,A q throughout most of the substrate; it is reasonable to 
assume that the A's have little influence, i.e., P«P T . Note that 
if the approximation works at all, it is not limited to locations 
where P is large. The right sides of (4.11) can be thought of as 
driving terms, analogous to charge density, which have accumulat- 
ing effects in the sense that the solution anywhere is influenced 
by the charge density everywhere. If the charge densities are 
nsarly equal throughout most of the substrate, where they are 
greatest, the solutions will be nearly equal everywhere. If 
P>>p o /2 throughout most of the substrate, so that P~P T throughout 
most of the substrate, we will also have P~P T near S 1 where P is 

small. 
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Quantitative calculations of P T can be done by using the spe- 
cial function F (discussed in Appendix C) to write (4.10) as 

P T = ( Po /2) F ( 2A/p c , 2(H+0)/p o ) (4.12) 


but approximations are useful for visualization. Note that (4.12) 
and (3.11) are the same except that fl is replaced by n+0. The 
generalized ambipolar approximation is obtained by making the 
same replacement. Neglecting (p 0 /2Vip)V 2 compared to P 2 for the 
negative V 2 case, the approximation is 

P T » n + 0 if V 2 < 0 and P 2 » p 0 /2. (4.13) 

If V 2 >0, there is an AR and HRR separated by an ARB, which is the 
constant n+0 surface characterized by 


n + 0 = (p 0 /2V T ) V 2 defines ARB . (4.14) 


The approximation in the AR is 


P T a n + 0 - (P q / 2V t ) V 2 in AR if V 2 >0 and P 2 »p Q /2 . (4.15) 

Quantitative estimates of P in the HRR (and anywhere else) should 
use (4.12); but, for visualization purposes, it is enough to know 
that P T is much smaller in the HRR than in the AR. 

Returning to I e 1# the approximation is obtained by replacing r 
in (4.7) with F T defined by 


div [ (P T + A) grad 

r T ] = 0 in substrate 

(4 . 16a) 

r v = 

0 on 

(4.16b) 

r T = 

1 on S 2 . 

(4.16c) 
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With surface integrals expressed in terras of R Q and G 1 , and A Q 
related to D*, the approximation can be written as 


I 6/1 « - ( D e / D h^ ( A /p 0 ) ^ V 2 + 2V T p 2/ p o^/ R o " G 1 


+ 2D e q [(A 0 - A)/p 0 ] 


sub 


grad r Y • grad <p d 3 x . (4.17a) 


Currents at S 2 are estimated by substituting the above result 
into (4.4) to get 


l hf2 * " (1 “ A/p 0 ) (V 2 + 2V T p 2/Po>/ R o + G 2 

grad r Y • grad <p d 3 x (4.17b) 


- 2D h q [(A q - A)/p 0 ) 


sub 


J e,2 ~ < D e/ D h) (A/Po) < V 2 + 2V T p 2 /P 0 )/ R o " G 2 


- 2D e q ((A 0 - A) / P 0 ] 


sub 


grad r T • grad <t> d 3 x . (4.17c) 


The equations in (4.17) are approximations,, but the particular 
combination of equations given by 

<D h /D e >Ie ,2 - I h ,2 “ < V 2 + 2V T P 2 /Po>/ R o - d + V D e> G 2 

is exact. 
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4 . 5 A Mathematical Theorem 


The integral in (4.17) has an interpretation (as a weighted 
average of grad<£) , but is difficult to numerically evaluate in 
three dimensions. The objective of this and the remaining sec- 
tions is to make (4.17) computationally manageable. The first 
step towards this objective is to derive a theorem relating 
volume integrals to surface integrals. The identity derived here 
is a little more versatile, for our applications, than the usual 
divergence theorem. 

Let S(v) denote the constant fi u surface characterized by n u =v. 
Note that v can be used as one coordinate in a curvilinear coor- 
dinate system. The value of v determines which constant n u sur- 
face a given space point lies on. Let r 3 and r 2 be two surface 
coordinates selected so that (T 1 ,r 2/ v) form an orthogonal system. 
If J is a sectionally continuous, but otherwise arbitrary vector 
field, we have 


J • grad n d 3 x 

sub 


J • grad n u h 3 h 2 h 3 dr 3 dT 2 dv 


where h 3 , h 2 , and h 3 are the scale factors for the coordinates 
r 1 , t 2 , and v, respectively. But h 3 is given by [4] 


h 3 = | grad flj" 1 


so the equation becomes 


J • grad n d 3 x 

J sub 



J • 


n h 3 h 2 


dr. 



dv 


where n is the unit vector in the direction of gradfl u . The double 
integral inside of the brackets is a surface integral on the n u =v 
surface, so the equation now becomes 
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j • ds dv (arbitrary J) . (4.19a) 


sub 


grad fl u d J x = 


S(v) 


The normal unit vector in the surface integral is in the direc- 
tion of increasing n u , so it is directed outward from the region 
between and S(v). We therefore have 


p 

% 

A p 

• 

Of 

W 

II 

1 

J • ds , 

j • ds — 

JS(0) J 

Si J 

S(l) 


J * ds 


A trivial generalization of the above steps gives 


|r(v) 


j • grad fl u d 3 x = 


v 

0 


j • ds dv 1 


(4.19b) 


S(v') 


where R(v) is the region between S x and S(v) 


4.6 A Special Family of Generation Rate Functions 

The second step towards the goal of making (4.17) computation- 
ally manageable is to confine our attention to a special family 
of generation rate functions. It will be assumed that g can be 
expressed in the form 


g = a(n u ) grad n u ♦ grad n u 


(4.20) 


for some function a. It is always possible to express g in the 
form (4.20) in one dimension because the product of the gradients 
is a constant and the argument of a is a linear function of the 
spatial coordinate. If the substrate has length L and we are 
given a g(x) with the origin selected so that S x is at x=0 and S 2 
is at x=L, then n u =x/L and a ( v) =L 2 g (vL) . But (4.20) imposes a 
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restriction in three dimensions. If we are selecting a g to 
represent a hypothetical case of our own choice, we can always 
select it to have the form (4.20). A more probable situation is 
one in which a g has been given and there is no a satisfying 
(4.20). We then look for an a that gives some kind of best fit, 
or at least a good fit (if possible). It is left to the user to 
find a fitting function a, but some guidance is given below. 

Selection of a fitting function a may be a little easier if a 
is related to familiar physical quantities. Such quantities are 
Gf, G 2 , and the volume integral of g. We start with 


<p = (n u /D*) 

*i 

*V 

a(v,) dv, dv - (1/D*) 

’ n u ’ 

* 

0 

0 

0 


dv^ dv 


(4.21) 


which can be verified by substituting (4.20) and (4.21) into 
(2.20). The gradient is given by 


D* grad 0 = 



”1 

'v 

’ n u 


0 

afv-iJ dv 1 dv - 

0 

a(v) dv 
0 






grad n u 


(4.22) 


so 


G 2 


q d’ 


grad <p • ds 




'1 

'1 

*v 

V T (°h Po V*" 1 


a(v) dv - 


a(v 1 ) dv x dv 


• 

0 

0 

0 


(4.23a) 


where we have used (2.21). Similarly, 


G, = V, 


T (°h Po R o) 


-1 


0 


afv,) dv, dv 


(4.23b) 


43 



and 


g d 3 x = G 1 + G 2 = (D h P D R 0 ) 


sub 


-1 


a(v) dv . (4.23c) 


The three equations (4.23) relate a to familiar physical quan- 
tities and may provide some guidance for those looking for a fit- 
ting function a (one good method is derived in Section 6.4). But 
the analysis given here goes in the other direction. It is as- 
sumed that a has been provided and the objective is to calculate 
other quantities from it. When going in this direction, it is 
convenient to express quantities in terms of B instead of a, 
where B is defined by 


B (v) = v 

*1 

" v 2 

a^) dv 1 dv 2 - 

r% 

V 

% 

0 

0 J 

0 


a(v 1 ) dvj^ dv 2 (4.24a) 


so that 


D* (p 0 /2) B’(v) 


’1 f v 2 

afv^ dv x dv 2 

J 0 J 0 


‘v 

a(v 1 ) dVj^ . 

J 0 


(4.24b) 


The only thing that we need <x for is to construct B and B' . The 
latter functions will be used from now on. Combining (4.24) with 
the previous equations gives 


0 = (P 0 /2) 

fi(n u ) 

(4.25) 

= [D e /(D e + D h )] 

(V t /R q ) B'(0) 

(4.26a) 

= - [D e /( D e + D h)3 

(V t /R q ) B'(l) . 

(4.26b) 


Another important quantity is the sum n +<p which is expressed as 
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n + <P - (P 0 /2) %(n u ) 


(4.27) 


where B m is a modified 6 defined by 


B m ( v ) = B(v) + (V 2 /V T + 2P 2 /p 0 ) V (4.28) 

and is trivially related to 6. A separate symbol is used only for 
notational brevity. We can write (4.12) in terms of B m as 

P T = (p 0 / 2 ) F( 2 A/ Po , B m (n u )) . 


For notational brevity, we will leave out the first argument and 
write the equation as 


P T = (p 0 /2) F(B m (n u )) (abbreviated notation) . (4.29) 


The integral in (4.17) can be evaluated by using (4.19) togeth- 
er with grad<^=(p 0 /2)B' (n u )gradn u to get 


grad r T • grad ep d 3 x = 

J sub 


(P 0 /2) 


’1 

B'(v) 

J 0 


grad r v 

. S(v) 


ds dv . 


(4.30) 


But 


U [(P 0 /2) F(B m (v)) + A]’ 1 dv 

J 0 

r T = (4.3i) 

‘1 

( (P 0 /2 ) F(B (v)) + A] -1 dv 

J 0 
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which can be verified by substituting (4.29) and (4.31) into 
(4.16). Taking the gradient of (4.31) and substituting it into 
(4.30) while using 

grad n u • ds = 

. S(v) 

gives 

grad • grad <p d^x = Vip (2q Dj^ Rq) (INT2/INT1) (4.32) 

sub 

with the two integrals INTI and INT2 defined by 

INTI ■ ((P 0 /2) F(6 m (v)) + A]" 1 dv (4.33a) 

J 0 

INT2 ■ f [ (p 0 /2) F ( ( v ) ) + A] -1 &’ (v) dv . (4.33b) 


The ratio INT2/INT1 is a weighted average of 6', similar to the 
weighted average of d0/dx in the one-dimensional equations (4.8) 

and (4.9). 

The currents are estimated by substituting (4.26) and (4.32) 
into (4.17) to get 

I h / 2 « - (1 " A/ Po ) (V 2 + 2 V t P 2 /P 0 ) / r o 

- (D e /D h ) (A 0 /p 0 ) (V T /R Q ) B'(D 

- t(A 0 - A)/ Po ] (V t /R q ) (INT2/INT1) (4.34a) 


grad n u • ds = V T (q D h p Q R Q ) 
• s 2 
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(4.34b) 


I e,2 * (°e/ D h) (A/P 0 ) ( V 2 + 2V T p 2 /Po)/ R o 

+ (D e /D h ) (A 0 /p 0 ) (V t /R 0 ) B'(l) 

- (D e /D h ) [ (A q - A)/ Po ] (V t /R q ) (INT2/INT1) 


(Dh/ D e> X e, 2 ~ p h,2 = ( v 2 + 2V T p 2 /P 0 )/ R o + (V T /R 0 )B*(1) . (4.35) 

The two equations (4.34) are approximations while (4.35) is 
exact. Any two of the above three equations can be used to solve 
for the currents. 


4 . 7 A Numerical Integration 

With a function 6 given, all quantities on the right sides of 
(4.34) can be calculated, but the integrals INTI and INT2 given 
by (4.33) require numerical methods. The numerical integration is 
regarded as part of the theory, rather than an exercise left for 
the reader, so some discussion is given here. 

The reader might notice that some of the integration can be 
done analytically. The derivatives 6* and B m ' differ by a con- 
stant, so both integrals can be evaluated if we can evaluate INTI 
and the integral 


[(Po/2) F(6m(v)) + A] -1 V(v) dv 


[ (P 0 /2) F + A]” 1 (dfi in /dF)dF . 


F(B m ) is related to B m by 


F + (1 - 2A/ Po ) ln(l + p Q F/2A) = B m 


which allows d6 m /dF to be expressed in terms of F alone. The 
above integral can be expressed in closed form, so only INTI 
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requires numerical approximations. This method is intentionally 
not used, because it is equivalent to approximating r T in (4.31) 
by retaining the numerator on the right while approximating the 
denominator with a numerical estimate. Any error in the estimate 
upsets the normalization condition r T (l)=l. The estimates of the 
currents are insensitive to errors in r or in r T when properly 
normalized, but estimates are sensitive to errors that disturb 
the normalization. If this method is used, accurate current 
estimates require an accurate numerical estimate of INTI. This is 
not easy, because the integrand can be extremely skewed, requir- 
ing a carefully selected variable step size for accurate numeri- 
cal approximation. It is desirable to eliminate the need for such 
numerical sophistication by using a different method to evaluate 
the integrals. 

One simple method is to numerically approximate both integrals, 
using the same step sizes for both. To see why this works, note 
that the ratio INT2/INT1 is a weighted average of 6'. Even if the 
weight factor is extremely skewed, the step size need be no 
larger than dictated by 6' (i.e., the step size only needs to be 
small enough for 6' to be nearly constant in each subinterval) if 
the numerical approximation of the weight factor is correspond- 
ingly skewed and normalized. By using the same step sizes for 
both integrals, we insure that the numerical approximation of the 
weight function is normalized, even if the step sizes are not 
small enough for an accurate estimate of INTI. We can therefore 
use a uniform step size to evaluate the integrals. 

One potential source of numerical error, which gets worse with 
smaller step sizes, can and should be avoided. This error source 
is the subtraction of nearly equal numbers that will occur when 
using B'dv=dB. It is better to leave B'dv as it is. This means 
that the user is required to supply B' in addition to B, but this 
is not a lot of extra work. If the user can calculate 6 from 
(4.24a), than the user can also calculate B' from (4.24b). 

A suggested numerical integration is the following. Select a 
moderately large value for M (the numerical examples in Chapter 6 
used M=100) and then calculate the quantities listed below 
(arrays are obviously unnecessary if the quantities are calculat- 
ed when needed) : 
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*1 = 2A/ Po 

B^ 1 = B 1 (i/M) i=0, . . ,M 

B l = B ( i/M) + (V 2 /V T + 2P 2 /p 0 ) (i/M) i=O f ...,M 

Ci = C(P 0 /2) F^B^ + A]" 1 i=0,...M 

INTI « (C Q + C m )/( 2M) + (1/M) “s 1 C ± 

1=1 ■ L 

INT2 a (C Q B 0 ' + C M B M ')/(2M) + (1/M) C ± B^ 
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5. THE COMPLETE SOLUTION 


5 . 1 Introduction 


This chapter does little more than list the equations in Chap- 
ter 4 together with those in Appendix A, to produce a complete 
equation set that is able to solve for all currents and boundary 
values. The only effort required here is associated with nuisance 
details such as including an electrode-semiconductor contact 
potential, and selecting a notation common to both equation sets. 
Following the list of equations is a suggested algorithm for 
constructing device I-V curves. This algorithm is interpreted as 
the "complete solution." A simple necessary condition for satura- 
tion is derived in the last section. 


5.2 Notation 


The notation used for the substrate analysis is familiar by now 
and the notation used for the DR analysis is listed in Appendix 
A. Redundant notations are related below so that the redundancy 
can be eliminated. The scalar current densities in the DR equa- 
tions are evaluated at the DRB on the lightly doped side, which 
is S 2 . These currents are positive when directed from the n-side 
towards the p-side, so 


I h,2 “ Jh a D' T e , 2 
I h,2 = ^h A D' I e,2 


- j e A d for the p-type substrate 
j e A d for the n-type substrate 


where A^ is the DRB surface area. The total current I is also 
taken to be positive when directed from the n-side towards the p- 
side, so 


I T,2 ” 


I T,2 


I = - j T a d for the p-type substrate 
I = j T A d for the n-type substrate 

j T = 3 h + j e 


(5.1a) 

(5.1b) 

(5.1c) 
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The equilibrium majority carrier density is equated to the doping 
density, so 


p 0 = N a for the p-type substrate 

n o = N D for the n_t yP e substrate. 


The equilibrium minority carrier density was left out of the 
substrate equations, but retained in some of the DR equations. We 
therefore use 


jip as p 2 + n o where n Q = n i 2 /N A for the p-type substrate 

p n = p 2 + P 0 where p Q = n i 2 /N D for the n-type substrate 


where n^ is the intrinsic electron density. 

Contact potentials between electrodes and semiconductor are 
simulated by fictitious power supplies of voltage V c as shown in 
Figure 5.1. The p- and n-type substrates are both shown. In each 
case, the polarity of the fictitious power supply is chosen so 
that V c is positive. V c is given by the well-known equation 

V c = V T ln(N A Njj/ni 2 ) . (5.2) 

Lumped resistors R c (Figure 5.1) simulate ohmic contact resist- 
ances, and may also include any other desired circuit resistances 
associated with electrical connections outside of the diode 
interior. The voltage V is applied to the upper contact (Figure 
5.1) above any resistor elements that are included in R c . A 
current arrow indicates the direction of the current when I is 
positive, consistent with the sign convention stated above. The 
potentials are related by 


52 


V 


V 



Figure 5.1: Qualitative sketch of both diode types showing R c and 
Vq. The currents are positive when in the indicated directions. 
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V 2 + V DR - V C + 1 R C 

for 

p-type 

substrate 

(5.3a) 

v 2 “ V DR + V C " 1 R C 

for 

n-type 

substrate. 

(5.3b) 


5.3 Equation Summary for the n + /p Diode 

All equations, excluding those listed in Section 5.2 and geo- 
metric information that must be supplied by the reader, are 
listed here for the p-type substrate diode. 

Starting with the doping densities N A (p-side) and Nq (n-side) , 
the low field mobilities M 0> h and ^o,e' the saturat i° n velocity 
v, the thermal voltage V T , the elementary charge q, and the di- 
electric constant e, other constants are calculated from 


D h =V T M 0 , h . D e - V T Mo . e , D* = 2D h D e /(D h + D e ) (5.4a) 


'e ~ V T **o,e 


o Q - (q/V T ) D h N A , A c = D h N A / (D h + D e ) 


(5.4b) 


a e = 1/ (q V T M 0 , e > ' v T b = V & v > * 


(5.4c) 


Boundary values P 2 end V 2 must be solved. The parameters A and E 
are defined in terms of P 2 and V 2 by 


E = [ln(l + 2P 2 /N a )] _1 if V 2 = 0 and P 2 > 0 . (5.5a) 


If v 2 +0 and V 2 + -2 V t P 2 /N a , use 


Z± = (V T /V 2 ) (1 + 2P 2 /N a ) , Z 2 ^ V T /V 2 (5.5b) 


E = H(Z 1 ,Z 2 ) 


(5.5c) 


where the special function H is defined in Appendix B. For any 
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case such that V 2 =f -2V T P 2 /N A , use either 


A = (N a /2) [1 - (V 2 /V T ) E] 


( 5 . 5d) 


or 


A = P 2 e 1 / E (1 - e" 1 ^) -1 


( 5 . 5e) 


The two equations for A give the same result in theory, but the 
second should be used if (V 2 /V T )E is so nearly equal to 1 that 
the first requires more numerical precision than is available. 
Otherwise, the first can be used. The functions n u , 0, and n are 
defined by 


div grad n u =0 in sub. , n u =0 on electrode, fl u =l on DRB 
div grad 0 =-g/D* in sub. , 0=0 on electrode, 0=0 on DRB 


n = [P 2 + (n a / 2 V t ) v 2 ] n u 


with reflective boundary conditions on the insulated boundaries 
tacitly assumed. The electrical resistance between electrode and 
DRB produced by the uniform conductivity a Q is R Q . The ambipolar 
diffusion currents and G 2 are given by 


Gi = - q D 


grad 0 * ds 


(i = 1,2) 


with the unit normal vector chosen so that G^ is positive. R Q and 
the G's may depend on the DR width W. 

An approximation for P applicable when P»N A /2 throughout most 
of the substrate is P T given by 
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P T = (N a / 2) F(2A/N a , 2 (n+0) /N a ) 


where the special function F is discussed in Appendix C. The 
simpler generalized ambipolar approximation is useful for visual- 
ization when P 2 »N a / 2. If V 2 <0, the approximation is 

P T « n + 0 if V 2 < 0 and P 2 » N A /2. 

If V 2 >0, there is an AR and HRR separated by an ARB, which is the 
constant n+0 surface characterized by 

n + <p = (N a / 2V t ) V 2 defines ARB . (5.6) 

P T is small in the HRR, but the approximation in the AR is 

P v « n + 0 - (N a / 2V t ) V 2 in AR if V 2 >0 and P 2 »N A /2 . 

Approximations for the currents are obtained by first defining 
r T by 

div ( (P T + A) grad r T ] =0 in substrate 
r T = 0 on electrode, r T =1 on DRB . 

The currents are approximated by 

j h A d « (1 - A/N a ) (V 2 + 2V T P 2 /N A )/R 0 - G 2 

+ 2D h q C (A 0 - A)/N a ] grad r T • grad <p d 3 x 

J sub 

(D h /D e )j e = j h +(1+ D h /D e )G 2 /A D - (V 2 + 2V T P 2 /N a )/(A d R q ) . (5.7) 
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Calculations are manageable in three dimensions if g can be 
expressed as 


g = o(n u ) grad n u • grad n u 


(5.8) 


for some function a, which is used to construct the user— supplied 
function 6 and derivative 6' given by 


D (N a /2) 6 (v) = v 


*1 

V 2 

a(v-,) dv, dv 2 - 

v r 

0 

0 

0 


“(Vj) dV 2 dv 2 (5.9a) 


D (N a /2) B'(v) = 


a(v 1 ) dv 1 dv 2 - 


a(v 1 ) dv 1 . (5.9b) 


The modified B is given by 


B m (v) = B(v) + (V 2 /V T + 2P 2 /N A ) v (5.9c) 

so that 

n + <t> = (N A /2) B m (n u ) (5.10) 

P T = (N a / 2) F(2A/N a , B m (n u ) ) . (5.11) 

The currents are now approximated by 


j h * C (N A - A)/N a ] (V 2 + 2 V t P 2 /N a )/(A d R q ) 

+ V T [(N a - A 0 )/N a ] B'(1 )/(A d R q ) 

+ V T [(A 0 - A)/N a 1 (INT2/INT1 )/(A d R q ) (5.12a) 
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(D h /O e ) j e " ^h ” [V 2 + V T B'(l) + 2V T P 2 /N A ]/(A D R d ) (5.12b) 


where the two integrals INTI and INT2 are evaluated by selecting 
a moderately large M (e.g., 100) and using 



X x = 2 A/N a 

(5.13a) 


Bi * = B ' (i/M) i=0, . . ,M 

(5.13b) 

Bi = & ( i/M) 

+ (V 2 /V T + 2P 2 /N a ) (i/M) i=0 , . . . ,M 

(5.13c) 

c i “ 

[(N a / 2) F(X 1 ,B i ) + A]" 1 i=0,...M 

(5.13d) 


M-l 

(5 . 13e) 

INTI « (C 0 + C m )/( 2M) + (1/M) Ci 


M-l 

. ( 5 . 13f ) 

H 

Z 

to 

n 

o 

o 

B q 1 + C M B m ')/( 2M) + (1/M) Ci Bi' 

that (5.12b) 

can be rewritten as 


< 

to 

II 

> 

a 

0 

(j h " D h ^e/ D e) ” 2V T P 2/ N A " V T B ' (1) 

. (5.14) 


One of the DR equations is 


j h = q W g D 


(5.15) 


with W the DR width and g D the value of g at the DR location. 
Another DR equation is 
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exp(- Vqj^/ V tj.) — Nq [P 2 + n Q - Vijb jip] 

— Nq ^ (Vrp e/q) (a e jp) 2 [N A 
Vrpb 3 <p 3 ^ “ Vrpb 3 <p 2P 2 2n 0 ] ^ xf 3p ^ 0 

exp(- V DR /V T ) = N d -1 (P 2 + n Q ) if j T < 0 

which can be solved for P 2 in terms of j T and the DR voltage drop 
V DR using 


P 2 = N d exp(- V DR /V T ) - n Q if j T < 0 (5.16a) 

Tp = exp (- + Vrpb jrp (5.16b) 

T 2 = (V t / 2) (e/q) (a e j T ) 2 /(N A + V T b j T ) (5.16c) 

T 3 = N d exp (- V DR /V T ) + (1/2) V T b j T + N a /2 (5.16d) 

T 4 = T 3 + [T 3 2 + 4T 2 ] 1 / 2 (5 . 16e) 

P 2 = T i - n D + 2T 2 /T 4 if j T > 0 . ( 5 . 16f ) 


The DR equation used to solve for W is 

W = (2e/q) 1 / 2 VpR 1 / 2 [(N a + V T b j T ) V6 / 2 

+ (2e/q) 1 / V6 (Vrp a e j T ) V6/3 V DR -1 / V6 ] _1 / V6 if j T > 0 (5.17a) 

W = [ (2e/q) V dr /N a ] 1 / 2 if j T < 0 (5.17b) 
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5.4 Algorithm for Constructing the n + /p Diode I-V Curve 

A suggested algorithm for constructing I-V curves for the p- 
type substrate diode is listed below. The voltage polarity and 
direction of current when positive are shown in Figure 5.1. An 
example of an I-V curve is seen by looking ahead to Figure 6.1 in 
Chapter 6. The diode delivers power (solar cell operation) when V 
is negative (a forward-biasing polarity) with I positive (a 
reverse current produced by photogeneration) . The "model curve 
for the particular example shown in the figure saturates for V 
greater than about -0.4 volts. Numerical problems will result if 
yg ■try to extend the curve too far into saturation , because A 
calculated from (5.5) becomes so close to zero that finite numer 
ical precision fails to distinguish it from zero. But there is no 
need to extend the plot beyond the point where such a problem 
first occurs, because such a point is far into saturation. In the 
opposite extreme of small (negative) V, the curve is very steep. 
Attempting to extend the curve too far in this direction also 
produces numerical problems because some calculated guantities 
become extremely sensitive to tiny errors (smaller than machine 
precision) in other guantities. But there is no need to extend 
the plot beyond the point where such problems begin to occur, 
because the current is large enough (in absolute value) to de- 
stroy the device. The objective is to plot points in the "range 
of interest," which is the range that avoids numerical problems 
and should also be the range that is physically interesting. A 
suggested algorithm is the following: 

2 

(1) Assign values to q, e/q, V T , N A , N D , n Q (=n^ /N A ) , V c 
(using (5.2)), R c , A D , g D , and the constants on the left 

sides of (5.4). 

(2) Select a positive value for Vp R . Each selected value 
will produce one point on the I-V curve. Trial and error is 
the simplest way to find a V DR value that produces a point 
in the range of interest. After several I-V points have been 
plotted, they can guide later selections of V DR values. 

(3) Guess at a value for j T . 

(4) Use (5.16) to solve for P 2 • Change the value to zero if 
the presence of n^ in (5.16) produces a negative value. 

(5) Use (5.17) to solve for W and (5.15) to solve for j^. 
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Then calculate j e from j T -j h . 


(6) With a value assigned to the DR width W, the substrate 
geometry is also specified. Assign a value to R Q . Find a 
fitting function a that (approximately) satisfies (5.8), and 
use (5.9) to construct the functions B and B'. 

(7) Use (5.14) to solve for V 2 . 

(8) Use (5.5) to solve for E. The function subprogram in 
Appendix B can be appended to any FORTRAN driver code for 
numerical evaluation of the function H. Note that the com- 
puter version of H contains a redundant argument Z 3 for 
improved numerical accuracy. Before calculating E, first 
calculate Z 3 from Z 3 =2 (V T /V 2 ) (P 2 / N A ) • Then calculate E from 
E=H (Z 1 ,Z 2 ,Z 3 ) . 

(9) Use (5.5) to solve for A. If A is found to be negative, 
the j T guess was probably too large. Try a less positive or 
a more negative j T . If A is positive but so close to zero 
that the available numerical precision cannot distinguish it 
from zero when (5.5e) is used, it is probable that either 
the j T guess was too small, or the V DR selection places the 
I-V point too far into saturation. First try a larger j T . If 
convergence (step 12 below) cannot be obtained with j T large 
enough to avoid this problem, use a smaller V DR . 

(10) Use (5.13) to calculate the integrals INTI and INT2 . 
The function subprogram in Appendix C can be appended to any 
FORTRAN driver code for numerical evaluation of the function 
F. 

(11) Use (5.12a) to calculate a new value for j h , denoted 
j h>new * Then calculate S j h s j hf n ew"V Calculate I from 

and then use (5.3) to calculate V. 

(12) Repeat steps 3 through 11 using different jip guesses 
until sufficiently close bracketing guesses have been found. 
Two guesses bracket the actual value if they produce Sj^'s 
having opposite signs. Bracketing guesses are sufficiently 
close when V and I calculated from the two guesses both 
agree, within the required precision. It is often necessary 
for bracketing guesses to have four- or five-digit agreement 
in order for the two V estimates to have three-digit agree- 
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ment. When the required agreement has been obtained, plot the 
I-V point and go back to step 2 for additional points. 


5.5 Equation Summary for the p + /n Diode 

All equations, excluding those listed in Section 5.2 and geo- 
metric information that must be supplied by the reader, are 
listed here for the n-type substrate diode. 

Starting with the doping densities (p-side) and Nq (n-side) , 
the low field mobilities h and M 0 ,e' the saturation velocity 
v, the thermal voltage V T , the elementary charge q, and the di- 
electric constant e , other constants are calculated from 


D h ~ v t ^o,h ' D e - V T ^o,e ' D - 2D h D e/ ( D h + D e) (5.18a) 
a 0 = (q/V T ) D e N D , A 0 = D e N D / (D h + D e ) (5.18b) 

a h = 1/ (q V T ju 0 ,h) ' v T b = */(q v > ‘ (5.18c) 

Boundary values P 2 and V 2 must be solved. The parameters A and E 
are defined in terms of P 2 and V 2 by 

E = [ln(l + 2P 2 /N d )] -1 if V 2 = 0 and P 2 > 0 . (5.19a) 

If v 2 =J=0 and V 2 =f 2 V t P 2 /N d , use 

Z 1 s " ( v t/ V 2 ^ + 2P 2/ N d) ' z 2 55 ~ v t/ V 2 (5.19b) 

E = H(Z 1 ,Z 2 ) (5.19c) 


where the special function H is defined in Appendix B. For any 
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case such that V 2 =j=2V T P 2 /N D , use either 


A = (N d / 2) [1 + (V 2 /V T ) E] ( 5 . 19d) 

or 

A = P 2 e _1 / E (1 - e -1 / 15 ) -1 . (5 . 19e) 


The two equations for A give the same result in theory, but the 
second should be used if (V 2 /V T )E is so nearly equal to -1 that 
the first requires more numerical precision than is available. 
Otherwise, the first can be used. The functions n u , 0, and D are 
defined by 


div grad n u =0 in sub. , n u =0 on electrode, 
div grad 0 =-g/D* in sub. , 0= 0 on electrode 

n = [P 2 - (n d /2v t ) v 2 ] n u 


n u =l on DRB 
0=0 on DRB 


with reflective boundary conditions on the insulated boundaries 
tacitly assumed. The electrical resistance between electrode and 
DRB produced by the uniform conductivity a Q is R Q . The ambipolar 
diffusion currents and G 2 are given by 


_ _ * 
G-l = - q D 


grad 0 • ds (i = 1,2) 

S. 


with the unit normal vector chosen so that G^ is positive. R Q and 
the G's may depend on the DR width W. 

An approximation for P applicable when P»N D /2 throughout most 
of the substrate is P T given by 
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P T = (N d /2) F ( 2A/N d , 2(n+0)/N D ) 


where the special function F is discussed in Appendix C. The 
simpler generalized ambipolar approximation is useful for visual 
ization when P 2 »N D /2. If V 2 >0, the approximation is 

P T « fl + <p if V 2 > 0 and P 2 » N D /2. 

If V 2 <0, there is an AR and HRR separated by an ARB, which is the 

constant n +<p surface characterized by 

n + <p = - (N d / 2V t ) V 2 defines ARB . (5.20) 

P T is small in the HRR, but the approximation in the AR is 

P T « fl + <P + (N d / 2V t ) v 2 in AR if V 2 <0 and P 2 »N D /2 . 

Approximations for the currents are obtained by first defining 
r T by 

div [(P T + A) grad r T ] = 0 in substrate 
r T = 0 on electrode, r T = 1 on DRB . 

The currents are approximated by 


j e A D « (1 - A/N d ) (2 V t P 2 /N d - V 2 )/R q - G 2 


+ 2D e q [ (A q - A) /N d ] 


grad F T 

J sub 


grad <p d 3 x 


(D e / D h) jh = ^e + (1 + D e/°h) G 2 /A D " (2V T P 2/ N D ” V 2)/( A D R o^ * 
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Calculations are manageable in three dimensions if 
expressed as 


g = <*( n u ) grad n u • grad n u 


for some function a, which is used to construct the user 
function B and derivative B* given by 


D*(N d /2) B (v) = v 

'1 

V 2 

a(v,) dv, dv, - 

*v r 

• 

0 

0 

0 


afv-j^ dv-^ dv 2 


D (N d /2) B'(v) = 


afv-j^) dv ^ dv 2 - 


a(v x ) dv-L 


The modified B is given by 


6 m( v ) = B(v) + (2P 2 /N d - V 2 /V T ) v 

so that 

n + * “ < n d/ 2) B m< n u> 

P T - (N D /2) F(2A/N d , B m (n u ) ) . 
The currents are now approximated by 


g can be 

(5.21) 

-supplied 

(5.22a) 

(5.22b) 

(5.22c) 

(5.23) 

(5.24) 
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(5.25a) 


j e « [ (N d - A)/N D ] (2V t P 2 /N d " v 2)/( A D R o> 

+ V T [ (N d - A q )/N d ] B'(1 )/(A d R q ) 

+ V T [ (A q - A)/N d ] (INT2/INT1 )/(A d R q ) 


<D e /D h ) 3h - U - [ 2V<j> P 2 /N d + V T 6MD - V 2 ]/(A D R 0> (5 . 25b) 

where the two integrals INTI and INT2 are evaluated by selecting 
a moderately large M (e.g., 100) and using 



X 1 = 2A/N d 

(5.26a) 


B^ 1 = G 1 (i/M) i=0 , • • f M 

(5.26b) 

B i 

= B (i/M) + (2P 2 /N d - V 2 /V T ) (i/M) i=0,...,M 

(5.26c) 


C i = [(N d / 2) F(X 1 ,B i ) + A] -1 i=0,...M 

( 5 . 26d) 


M-l 

INTI » (C 0 + C m )/( 2M) + (1/M) .S i C ± 

( 5 . 26e) 

INT2 * (C Q B 0 ' + C M B m ')/(2M) + (1/M) C ± B^ 

. ( 5 . 26f ) 

Note that 

(5.25b) can be rewritten as 


V 2 = 

- a d r o (j e ” D e 3h/ D h^ + 2V T p 2/ N D + V T fi ' (1) 

. (5.27) 

One of 

the DR equations is 



j e = <3 w % 

(5.28) 


with W the DR width and g D the value of g at the DR location 
Another DR equation is 
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exp(- V DR /V T ) = N a _1 [P 2 + P 0 - V T b j T ] 

- N a _1 (V t e/q) (a h j T ) 2 [Np 
+ V T b j T ] _1 [N d - V T b j T + 2P 2 + 2p 0 ] _1 if j T > 0 

exp(- V DR /V T ) = N a _1 (P 2 + Po ) if j T < 0 

which can be solved for P 2 in terms of j T and the DR voltage drop 


v dr usin< ? 

P 2 = N a exp (- V DR /V T ) - p Q if j T < 0 (5.29a) 

T 1 “ N A ex P(~ V DR/ V T) + V T b (5.29b) 

T 2 = (V t / 2) (e/q) (a h j T ) 2 /(N D + V T b j T ) (5.29c) 

T 3 = N a exp (- V DR /V T ) + (1/2) V T b j T + N d /2 (5.29d) 

T 4 = T 3 + [T 3 2 + 4T 2 ] 1 / 2 ( 5 . 29e) 

P 2 = T i ~ P 0 + 2T 2 / T 4 if j T > 0 . (5 . 29f ) 


The DR equation used to solve for W is 

W = (26/q) 1 / 2 V^ 1 / 2 ((N D + V T b j T ) V6/2 

+ (2e/q) 1 / V6 (V T a h j T ) V6 / 3 V DR " 1 / V6 ] " 1 / V6 if j T > 0 (5.30a) 

W = [ (2e/q) V^/Np] 1 / 2 if j T < 0 (5.30b) 
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5.6 Algorithm for Constructing the p + /n Diode I-V Curve 

A suggested algorithm for constructing I-V curves for the n- 
type substrate diode is listed below. The voltage polarity and 
direction of current when positive are shown in Figure 5.1. An 
example of an I— V curve is seen by looking ahead to Figure 6.8 in 
Chapter 6. The diode delivers power (solar cell operation) when V 
is positive (a forward-biasing polarity) with I positive (a 
reverse current produced by photogeneration). The "at 3.0 /xm” 
curve for the particular example shown in the figure saturates 
for V less than about -0.3 volts. Numerical problems will result 
if we try to extend the curve too far into saturation, because A 
calculated from (5.19) becomes so close to zero that finite 
numerical precision fails to distinguish it from zero. But there 
is no need to extend the plot beyond the point where such a 
problem first occurs, because such a point is far into satura- 
tion. In the opposite extreme of large V, the curve is very 
steep. Attempting to extend the curve too far in this direction 
also produces numerical problems because some calculated quanti- 
ties become extremely sensitive to tiny errors (smaller than 
machine precision) in other quantities. But there is no need to 
extend the plot beyond the point where such problems begin to 
occur, because the current is large enough (in absolute value) to 
destroy the device. The objective is to plot points in the "range 
of interest,” which is the range that avoids numerical problems 
and should also be the range that is physically interesting. A 
suggested algorithm is the following: 

(1) Assign values to q, e/q, V T , N A , N D , p Q (=n i 2 3 4 /N D ) , V c 
(using (5.2)), R c , A D , g D , and the constants on the left 
sides of (5.18). 

(2) Select a positive value for Vp R . Each selected value 
will produce one point on the I-V curve. Trial and error is 
the simplest way to find a V DR value that produces a point 
in the range of interest. After several I-V points have been 
plotted, they can guide later selections of V^ R values. 

(3) Guess at a value for j T . 

(4) Use (5.29) to solve for P 2 . Change the value to zero if 
the presence of p Q in (5.29) produces a negative value. 
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(5) Use (5.30) to solve for W and (5.28) to solve for j e . 
Then calculate from j T ~j e . 

(6) With a value assigned to the DR width W, the substrate 
geometry is also specified. Assign a value to R Q . Find a 
fitting function a that (approximately) satisfies (5.21), 
and use (5.22) to construct the functions R and B'. 

(7) Use (5.27) to solve for V 2 . 

(8) Use (5.19) to solve for E. The function subprogram in 
Appendix B can be appended to any FORTRAN driver code for 
numerical evaluation of the function H. Note that the com- 
puter version of H contains a redundant argument Z 3 for 
improved numerical accuracy. Before calculating E, first 
calculate Z 3 from Z 3 =-2 (V!p/V 2 ) (P 2 /Nq) . Then calculate E from 
E=H (Z lf Z 2 ,Z 3 ) . 

(9) Use (5.19) to solve for A. If A is found to be negative, 
the j T guess was probably too large. Try a less positive or 
a more negative j T . If A is positive but so close to zero 
that the available numerical precision cannot distinguish it 
from zero when (5.19e) is used, it is probable that either 
the jip guess was too small, or the selection places the 
I-V point too far into saturation. First try a larger j T . If 
convergence (step 12 below) cannot be obtained with jip large 
enough to avoid this problem, use a smaller V^. 

(10) Use (5.26) to calculate the integrals INTI and INT2 . 
The function subprogram in Appendix C can be appended to any 
FORTRAN driver code for numerical evaluation of the function 
F. 


(11) Use (5.25a) to calculate a new value for j e , denoted 
j e , new * Then calculate 6 j e =j S/ new “j e - Calculate I from j.pA D 
and then use (5.3) to calculate V. 

(12) Repeat steps 3 through 11 using different j T guesses 
until sufficiently close bracketing guesses have been found. 
Two guesses bracket the actual value if they produce <Sj e 's 
having opposite signs. Bracketing guesses are sufficiently 
close when V and I calculated from the two guesses both 
agree, within the required precision. It is often necessary 
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for bracketing guesses to have four- or five-digit agreement 
in order for the two V estimates to have three-digit agree- 
ment. When the required agreement has been obtained, plot 
the I-V point and go back to step 2 for additional points. 


5.7 A Necessary Condition for Saturation 

"Saturation" is defined here to mean that the diode current is 
virtually the same as the total rate that charge is liberated in 
the device via photogeneration. Looking ahead to Figures 6.3 and 

6.8 in Sections 6.2 and 6.3, we see that some I-V curves display 
saturation while others do not. Now that the DR and substrate 
equations have been listed together, we can derive a very simple 
necessary (but not sufficient) condition for saturation. Satura- 
tion, strong funneling, a wide HRR, and DR collapse occur togeth- 
er, so the condition derived below can also be regarded as a 
necessary condition to collapse a DR. 

We start with the n + /p diode where saturation means 


A D 3 e 


q 


g d 3 x = G^ + G 2 

J sub 


(5.31) 


where we have used (2.22b) . Using (5.31) and the DR equation 
(5.15) with the substrate equation (5.7) gives 


G 2 + q A d W g D « (D h /D e ) G x + (V 2 + 2V T P 2 / n a)/ R o 


which is a necessary and sufficient condition for saturation, but 
contains unknown boundary values. The only additional information 
regarding the DR needed to obtain a simpler necessary condition 
is the fact that the quantity 


V 2 + 2V T P 2 /N a 
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is positive. This quantity is obviously positive if V 2 is posi- 
tive. If V 2 is negative, we have forward-biasing conditions and 
P 2 / n a will be much larger than -V 2 /V T . We may therefore assume 
that the quantity is positive and the necessary condition becomes 


G 2 + q A D W g D > ( D h/°e )G l (necessary to saturate n + /p) . (5.32a) 


The left side of (5.32a) is the rate carriers are generated in 
the DR plus the rate that carriers flow into the DR as predicted 
by the ambipolar diffusion equation with homogeneous boundary 
conditions. On the right side, G ± is the rate carriers flow to 
the electrode as predicted by the same equation. The necessary 
condition states that the rate carriers are generated in the DR 
or flow into the DR must exceed a certain multiple of the rate 
they flow to the electrode, as predicted by ambipolar diffusion. 
This is a statement regarding the spatial distribution of photo- 
generation and says nothing about the strength of the photogener- 
ation. The condition is satisfied if carrier generation is con- 
fined to locations sufficiently close to the M J . This is clearly 
not a sufficient condition because it can be satisfied under 
LILC. But if the condition is not satisfied, the DR will not 
collapse even if the generation rate is great enough to result in 
p 2» N A' ijnplyi 1 " 1 ? that the latter condition is not sufficient to 
collapse a DR. This assertion is supported by computer simulation 
results discussed in Section 6.3. 

The analog of (5.32a) for the p + /n diode is 


G 2 + q A D W g D > (Dg/Dj^G-L (necessary to saturate p + /n) . (5.32b) 

Because D e >D h , (5.32b) is more difficult to satisfy than (5.32a). 
DR collapse requires carrier generation to be closer to the MJ 
for the p + /n device than required for the n + /p device. This is 
our first indication that funneling is more difficult to induce 
in the p + /n device. But (5.32a) and (5.32b) are only necessary 
(not sufficient) conditions and we cannot yet rigorously conclude 
that the p + /n device is less susceptible to funneling, although 
it is, as will be seen in Chapter 6. 
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6. NUMERICAL EXAMPLES AND CONCLUSIONS 


6.1 Introduction 


This chapter presents numerical examples to illustrate concepts 
already discussed and to inspire additional discussion. Unneces- 
sary complexity does not help here, and the examples will be 
simple. Sections 6.2 and 6.3 treat one-dimensional n + /p and p + /n 
diodes. Section 6.4 treats a simple three-dimensional problem 
having rotational symmetry. Conclusions are summarized in Section 
6.5. Qualitative sketches in Figure 5.1 (Chapter 5) show the 
polarity convention and the direction of the current when posi- 
tive. The n + /p diode delivers power (solar cell operation) when V 
is negative (a forward-biasing polarity) with I positive (a 
reverse current produced by photogeneration) . The p + /n diode 
delivers power when V is positive (a forward-biasing polarity) 
with I positive (a reverse current) . Readers that are not inter- 
ested in mathematical theory can ignore the paragraphs in the 
sections below that discuss & and J3'. 

Comparisons are made between theoretical (or model) predictions 
and predictions from a computer simulation code called PISCES 
[5] . Material constants used for the calculations are either 
default values used by PISCES or are derived from such values. 
All examples below used the following data (see Sections 5.2, 
5.3, and 5.5 for notation): 

doping density (substrate side) = 8 x 10 14 /cm 3 
doping density (other side) = 10 20 /cm 3 

Rc = 0 

= 1.5 x 10 10 /cm 3 
V T = 0.016 V 

q = 1.6 x 10~ 19 C, e/q = 6.536 x 10 6 /V-cm 
D h = 13/0 cm 2 /s, D e = 26.0 cm 2 /s 

a h = 4.84 x 10 17 /A-cm 2 , a e = 2.42 x 10 17 /A-cm 2 

V T b = 3.7 x lO^/A-cm 
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PISCES includes a variety of second-order effects, such as band- 
gap narrowing, several types of recombination mechanisms, and 
mobilities that depend on a variety of things. Good agreement 
between model and PISCES predictions indicates that the second- 
order effects are not important to the quantities of interest in 
the particular example considered. 


6.2 The One-Dimensional n + /p Diode 

We start with the one-dimensional n + /p diode. Let L be the 
distance between the electrode and MJ, so L-W is the distance 
between electrode and DRB, where W is the DR width. Two types of 
generation rate functions are considered. One is uniform below 
the MJ, i.e., g=g Q where g Q is a constant. The total rate per 
device area that carriers are generated below the MJ for this 
case is g Q L. For the other case, all carrier generation is con- 
fined to a horizontal plane a specified distance x Q above the 
electrode, so g=g 0 L$ (x-x Q ) where S is the Dirac delta function, x 
is the distance from the electrode, and g Q L "the total rate per 
device area that carriers are generated below the MJ . 

The only quantities used in the model that depend on geometry 
and/or carrier generation are Ap, g^, R Q , and the functions B and 
B* . The DRB area A D is also the device area and is set equal to 
lcm 2 , so that the device current in amps is also the current 
density in amps/cm 2 . For all cases, we use 


R q = (L - W)/(A d a Q ) . 


For the uniform case, we have g D =g Q and 


B(v) = [ (L - W) 2 /(N a D*)] g Q (1 - v) v 


B ' (v) = t(L - W) 2 /(N a D*)] g Q (1 - 2v) . 
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For the delta function case (with generation below the DRB) , we 
have g D =0 and 


B(v) = [2 L/(N a D*)] g D (L - W - x Q ) 
B(v) = [2 L/(N a D*)] g c x Q (1 - v) 

B’(v) = [2 L/(N a D*)] g Q (L - W - x 
(v) = - [2L/(N A d*) 3 g 0 x 0 


v if v < x Q / (L - W) 
if v > x q /(L - W) 

D ) if v < x q /(L - W) 
if v > x q /(L - W) . 


The above information supplements step 6 of the algorithm in 
Section 5.4. All other steps are explicit and require no addi- 
tional explanation. 

The dimension L is arbitrarily set equal to 5 urn in the exam- 
ples below. (It could be made larger but must be less than a 
diffusion length, because recombination is neglected in the 
substrate.) Examples are only interesting if they show signifi- 
cant deviations from classical theory predictions (implying high- 
injection— level— conditions) , and the generation rate was chosen 
to be large enough to make this happen. For this particular 
diode, a uniform generation rate of g=g 0 =l . 25x10^ ^/cm^-sec suf- 
fices. Including the factor of q, the total charge generation 
rate per device area below the MJ is 1000 amps/cm^ , which is the 
device current when saturated. 

Figure 6.1 compares model, PISCES, and classical predictions of 
the I-V curve produced by a uniform generation rate of 
1. 25xl0^^/cm^-sec, and shows that the classical prediction is not 
very good for this case. The classical prediction uses the clas- 
sical law of the junction, which is (5.16a) but used for all jip 
and with V DR set equal to V+V c - The classical estimate of W is 
used in (5.15) and to determine the electrode to MJ distance. The 
classical estimate is (5.17b) but used for all jp and with V^ R 
set equal to V+V c . The minority carrier substrate current is 
calculated by neglecting the drift term and calculating the 
carrier density from the minority carrier diffusion equation. It 
could be argued that classical theory is not being given a fair 
chance, because the ambipolar diffusion equation may be more 
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appropriate than the minority carrier diffusion equation for 
calculating carrier density. It turns out that the agreement in 
Figure 6.1 would be improved if the ambipolar diffusion equation 
was used for this calculation. But this improvement is not very 
satisfying in view of the fact that the very same "fix" will 
worsen the agreement for the p + /n diode under high-injection- 
level conditions treated in the next section. The best agreement 
obtainable from classical theory for the latter case is produced 
by the minority carrier diffusion equation. For consistency, this 
equation is used for all classical theory predictions. 

The model and PISCES predictions in Figure 6.1 show that satu- 
ration (I«1000 amps/cm 2 ) is reached even at some negative volt- 
ages. Saturation is an indication that funneling is very strong, 
but a better indication is obtained by looking at conditions 
(carrier density and voltage drops) inside of the device. The I-V 
point at V=1 volt is characterized by the following model-pre- 
dicted parameters: 


DR width (W) = 0.384 nm 
substrate voltage drop (V 2 ) = 1.627 volts 
electron density at DRB (P 2 ) = 9.171 x lO^/cm 3 
A = 6.123 x 10 -12 /cm 3 


A model-predicted estimate of electron density is P T (given by 
(5.11)), which is plotted from the above data against distance 
from MJ in Figure 6.2. The PISCES prediction is also shown. The 
PISCES prediction places the DRB closer to the MJ than the model 
prediction. (The DRB and ARB locations shown in the figure are 
model predictions.) This is consistent with the fact that PISCES 
calculates a smaller DR voltage drop (V^) than the model, and is 
probably due to band-gap narrowing, which PISCES includes but the 
model does not. Fortunately, this does not seem to affect the I-V 
curve in Figure 6.1. A compensating correction in the equilibrium 
built-in potential Vq allows PISCES and the model to agree on the 
device voltage drop V and the substrate voltage drop V 2 , even 
when they disagree on V^. If we account for the shift in DRB 
location, the two curves in Figure 6.2 will agree very well. 

A wide HRR is clearly shown in Figure 6.2, implying strong 
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Figure 6.2: Comparison of electron density predictions for the 
n + /p diode with a uniform g = 1.25 x 10 25 cm -3 s" 1 . 
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funneling. The V 2 for this configuration is 1.62 (PISCES) or 1.63 
(model) volts, which also implies strong funneling. About two- 
tenths of a volt is across the AR, with the remainder across the 
HRR, consistent with the statement that most substrate resistance 
is in the HRR. 

The effect of carrier generation location is interesting. The 
Figure 6.1 model curve is replotted in Figure 6.3, together with 
an I— V curve produced when all carrier generation is confined to 
a horizontal plane 2.5 jum above the electrode (more than 1 ura 
below the unperturbed DRB for biasing voltages up to 0.5 volts). 
The total generation rate below the MJ is the same for both 
cases. The I-V curves are so nearly identical that they could not 
be distinguished if smooth curves were drawn. Discrete points are 
shown to emphasize that there really are two data sets here, they 
just happen to lie on the same curve. It should not be concluded 
that the model predicts the uniform and "at 2.5 /m" cases to be 
equivalent. We can see differences if we look inside of the 
device, e.g., the DR and substrate voltage drops are individually 
different even when they have the same sums. Furthermore, classi- 
cal theory predicts a slightly larger current for the "at 2.5 nm" 
case. Therefore, there should be some difference between the two 
curves, but the difference is too small to be seen in the figure. 

Saturation in the "at 2.5 jum" curve implies that strong funnel- 
ing is induced at a distance, i.e., by carriers generated outside 
of the DR. To get the funneling process started, carriers must 
first diffuse to the DR. Once there, the DR partially collapses 
and a substrate electric field is created. This field drives more 
minority carriers to the DR and the funneling process becomes 
self sustaining . Figure 6.3 also shows the case where all carrier 
generation is 1 jxm above the electrode. Classical theory predicts 
a comparatively weak current for this case, because most carriers 
diffuse to the electrode where they recombine. The model shows 
that funneling is now diminished and no longer strong enough to 
produce saturation, but still strong enough for the current to be 
much larger than predicted by classical theory. 

Before ending this section, it should be verified that the 
model, PISCES, and classical predictions all come together under 
low— inj ection— level conditions. Such conditions are produced in 
the diode considered here by decreasing the carrier generation 
rate by two orders of magnitude. Figure 6.4 compares the predic- 
tions for the uniform but reduced generation rate and verifies 
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Figure 6.3: Comparison of model-predicted I-V curves for the n + /p 
diode when carrier generation location is varied. One curve is 
produced by a uniform g = 1.25 x 10 25 cm -3 s _1 (same as Fig. 6.1). 
For the other two curves, all carriers are generated at the 
indicated distance above the electrode. The total generation rate 
below the MJ is the same for all cases. 
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Figure 6.4: Comparison of I-V curve predictions for the n + /p 

diode with a reduced uniform g = 1.25 x 10 23 cm -3 s” 1 . 


81 



that the predictions do come together. 


6 . 3 The One-Dimensional p + /n Diode 

We now consider the same problem treated in the last section, 
except that n-type and p-type are interchanged. Figure 6.5 com- 
pares model, PISCES, and classical predictions of the I-V curve 
produced by a uniform generation rate of 1.25xl0^^/cm^-sec. The 
most noticeable difference between Figures 6.5 and 6.1 is that 
the p + /n diode is not saturating and the classical prediction is 
fairly good (although the classical prediction would not be as 
good if the ambipolar diffusion equation replaced the minority 
carrier diffusion equation, as discussed in the last section) . 
Compared to the n + /p diode under the same conditions, tunneling 
is greatly reduced for the p + /n diode. 

A closer comparison can be seen if the n + /p and p + /n curves are 
plotted on the same axis by replacing V with the bias voltage V B , 
where Vg=V for the n + /p diode and Vg=— V for the p^/n diode. In 
either case, reverse currents are positive and a positive Vg is a 
reverse-biasing polarity. The plot is shown in Figure 6.6. 

Classical theory predicts the p + /n device to have the larger 
(more positive or less negative) current at small Vg, with the 
curves coming together at larger Vg. This is understandable 
because the classical current is the sum of a forward current 
associated with biasing and a reverse current associated with 
photogeneration. The minority carrier currents, associated with 
photogeneration, at the electrode and DRB add up to the total 
generation rate in the substrate; the way this rate is divided 
between the currents at the two locations depends upon the spa- 
tial distribution of photogeneration, but not on mobility. The 
reverse current associated with photogeneration does not depend 
on mobility (mobility divides out of the equations) . But the 
forward current is reduced by a reduced minority carrier mobili- 
ty, so the device having the smaller minority carrier mobility 
(the p + /n diode) will have the larger net reverse current, unless 
the forward currents are negligible so that the two devices have 
the same currents. This is the classical prediction. 

The model prediction in Figure 6.6 agrees with the classical 
prediction in that the p + /n device has the larger current at 
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Figure 6.6: Comparison of n"*"/p and p"*"/n diode I—V curves with a 
uniform g = 1.25 x 10 25 cm -3 s _1 . 
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small V B . But for larger V B , funneling becomes stronger in the 
n + /p device and now this diode has the larger current. Depending 
on bias voltage and whether carrier generation is sufficient to 
produce funneling in at least one device, either device can have 
the larger current. 

When looking at either I-V points or I-V curves associated with 
different conditions, we may see a gradual transition between 
nonsaturation and saturation, and the two cases may not look so 
profoundly different. The two model points at Vg=l volt in Figure 
6.6 are not really very different. A more profound difference is 
seen if we look inside of the device at the carrier density and 
voltage drops. The p + /n point is characterized by the following 
model predicted parameters: 


DR width (W) = 1.123 » m 
substrate voltage drop (V 2 ) = -0.108 volts 
hole density at DRB (P 2 ) = 3.682 x 10 15 /cm 3 
A = 3.760 x 10 13 /cm 3 


The above data were used to plot the hole density in Figure 6.7, 
which also shows the PISCES prediction. PISCES predicts V 2 to be 
-0.113 volts, which is nearly the same as the model prediction. 
The agreement between the model and PISCES predictions looks good 
in Figure 6.7. 

Comparing Figure 6.7 and a V 2 value of about -0.11 volts to 
Figure 6.2 and a V 2 value of about 1.63 volts, we can now see 
striking differences between the two cases. The DR is collapsed 
and the substrate voltage drop is large for the n + /p case. But 
for the p + /n case, the DR is wide and supports nearly all of the 
applied plus built-in potential, with only a small fraction of 
this potential across the substrate. The n + /p case shows a wide 
HRR. There is a theoretically predicted HRR for the p + /n case, 
but it is so narrow as to be almost nonexistent. Because this HRR 
is so narrow, the substrate voltage is across a highly conductive 
region. This high conductivity nearly compensates for the small- 
ness of V 2 , so that funneling is occurring in this nonsaturated 
p + /n diode and the current is almost as large as in the saturated 
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Figure 6.7: Comparison of hole density predictions for the p + /n 
diode with a uniform g = 1.25 x 10^ cm ^ s 




n + /p diode. If the strength of funneling is measured by the size 
of the currents, then funneling is not greatly different for the 
two cases. If the strength of funneling is measured by the size 
of the substrate voltage drop (which is the convention used when 
strong funneling is equated to saturation) , then funneling is 
greatly different for the two cases. 

With the exception of a region close to the electrode, the 
minority carrier density in Figure 6.7 greatly exceeds the doping 
density, even at the DRB. It is interesting (perhaps surprising) 
that this is not sufficient to collapse the DR. The fact that the 
DR has not collapsed (enough for the substrate voltage to be 
great enough to produce saturation) can be predicted from the 
fact that the necessary condition (5.32b) is not satisfied. The 
condition can be satisfied if carriers are generated closer to 
the MJ. If all generation is moved to a horizontal plane 3 /m 
above the electrode, the necessary condition will be satisfied at 
any point on the I-V curve where the DR width W exceeds 0.5 /im. 
Assuming the generation rate is great enough to satisfy all other 
necessary conditions (whatever they are) , we can expect to see 
saturation somewhere on the I-V curve. This is seen in the "at 
3.0 /xm" curve in Figure 6.8. A close look at this curve finds a 
small but rapid change in slope at V«-0.3 volts. It seems 
reasonable to call this point the onset of saturation. The DR 
width near this point is between 0.78 and 0.74 jum (depending on 
the exact location of the onset point) , so the necessary condi- 
tion (5.32b) is fairly close to (but not quite) a sufficient 
condition for this example. 

Figure 6.8 also shows the I-V curve produced when all carrier 
generation is 2.5 fim above the electrode. The difference between 
this and the uniform case is large enough to be visible in the 
figure, but still very small. The "at 2.5 /xm" curve does not 
saturate, even though the generation location is only 0.5 /xm away 
from that for the saturating "at 3.0 /urn" curve. At V=-l volt, the 
substrate voltage drop for the "at 2.5 jum" case is about -0.11 
volt (almost the same as the uniform case), compared to -0.44 
volt for the "at 3.0 /urn" case. The currents for the two "at" 
cases are almost the same. This is another illustration of the 
fact that the difference between nonsaturation and saturation is 
more profound if we look at substrate voltage drops instead of 
currents . 

The final noticeable difference between Figures 6.1 and 6.5 is 
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Figure 6.8: Comparison of model-predicted I-V curves for the p + /n 
diode when carrier generation location is varied. One curve is 
produced by a uniform g = 1.25 x 10 25 cm -3 s” 1 (same as Fig. 6. 5). 
For the other two curves, all carriers are generated at the 
indicated distance above the electrode. The total generation rate 
below the MJ is the same for all cases. 



that the model does not agree as well with PISCES in the latter 
figure. This might be explained in terms of sensitivity. The 
condition represented in Figure 6.5 is close to some kind of 
threshold, in the sense that the device is trying to saturate but 
cannot quite do so. An HRR is wide enough to influence the minor- 
ity carrier current at the electrode, but not wide enough to 
either produce saturation or to be clearly visible in Figure 6.7. 
The calculated minority carrier current at the electrode is 
sensitive to error in the minority carrier density near the 
electrode where the density is small. It was argued in Section 
4.4 that if P t «P>>N d almost everywhere, than P T «P everywhere, 
even near the electrode. This is still true, but we must distin- 
guish P governed by the quasi-neutral equations from the PISCES- 
calculated minority carrier density, which is governed by a more 
complicated set of equations. While the model— and PISCES-pre- 
dicted minority carrier densities agree well in terms of absolute 
error, the relative or fractional error is significant near the 
electrode. We should expect some error when an HRR strongly 
influences the minority carrier current but does not block it, 
i.e., when conditions are almost but not quite able to produce 
saturation. But even under these adverse conditions, the agree- 
ment between the model and PISCES curves in Figure 6.5 is fairly 
good. 


6.4 A Simple Three-Dimensional Diode 

A simple three-dimensional example is considered, primarily to 
illustrate a general method for treating such problems. The 
objective is to illustrate the method while avoiding difficult 
integrals, so the example is highly idealized. Readers that are 
willing to evaluate difficult integrals can apply the method to 
more difficult problems. 

In this example, one DRB is isolated from all other DRBs . The 
DRB is a circular disk of radius r D and photogeneration is con- 
fined to a circular cylinder having the same radius r^ and length 
L. The cylinder is normal to the device and centered on the DRB. 
It is assumed that r^ and L are both small compared to the 
DRB-to-electrode distance. Because recombination is neglected, r D 
and L are both required to be small compared to the diffusion 
length. As long as the above conditions are satisfied, it is not 
required that the DRB-to— electrode distance be small compared to 
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the diffusion length. We can neglect recombination and regard the 
electrode as infinitely far away, so the problem to be solved 
reduces to that shown in Figure 6.9, which also shows the coordi- 
nate system. The generation rate is uniform and equal to g Q (a 
constant) inside of the cylinder. Cases in which the cylinder 
radius is less than r D might be approximated by the case consid- 
ered here if g Q is selected to produce the same total generation 
rate per unit length in the vertical direction. The special 
choice of r D for the cylinder radius simplifies some integra- 
tions. A better representation of a possible physical arrangement 
would use a generation function that is exponentially attenuated 
in the vertical coordinate. The attenuated problem is left for 
any reader that is willing to evaluate the required integrals. 

The DR width W is simulated by retaining the flat disk geometry 
but reducing the generation cylinder length from L to L— W (assum- 
ing that L>W) . A majority carrier current calculated from (5.15) 
or (5.28) compensates for the missing cylinder section. For 
notational brevity, a length L is used in the analysis and then 
replaced with L-W in the final equations. The DR width can also 
add to the lateral dimension r D in three dimensions, but this is 
ignored in the analysis below. No distinction is made between the 
DRB radius and the MJ radius. 

The only quantities, used in the algorithms in Sections 5.4 and 
5.6, that depend on geometry and/or carrier generation are A D , 
g D , R q , and the functions 6 and B'. We obviously have g 0 =g o and 
A D =irr D 2 . R 0 is well known for the flat circular disk and given by 
r =1/ (4a Q r D ) . The analysis is finished when the functions B and 
B* have been constructed. These functions are derived from a 
satisfying (5.8). But there is no such a for this three-dimen- 
sional problem and fitting is required. The definition of a "best 
fit" is somewhat arbitrary, but a particular definition will 
produce exact calculations of the ambipolar diffusion currents 
and G 2 * This is demonstrated below for arbitrary geometries and 
generation functions. Readers that are not interested in mathe- 
matical theory can go directly to the paragraph following the 
equations for B and B' on page 95. 

A sufficient condition for a fit to g to produce the correct 
G's is found by using (2.19), (2.20), and the divergence theorem 

to write (2.22a) as 
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Figure 6.9: A simple three-dimensional geometry. Generation is 
confined to the cylinder of length L and radius r^ (same as the 
DR radius) . The electrode is at infinity. The fl u = v surface en- 
closes the region R c (v) . 
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-Go/(q D ) = 


grad <p • ds = 


n u grad <p ♦ ds + 


n u grad <f> * ds 


<t> grad n u • ds - 


0 grad fl u • ds 


sub 


n u div grad <p d 3 x 


or 


Go = q 


sub 


n u g d 3 x 


We can select a number M and partition the substrate into subre- 
gions SjR, 6 2 R, . .., 5 m R, where each <SjR is the region between 
the n u =(i-l)/M and the n u =i/M surfaces. The integral can be 
written as the sum 


g d 3 x 


G 2 = 


M 

q 2 

i=l 


n 


5,R 


u 


g d 3 x 


M 

2 (i/M) 

i=l 


with the approximation on the far right becoming exact in the 
large M limit. A fitting function will produce the same G 2 
if it satisfies 


» 

c 

. S ± R 


gf it d3x = 


g d 3 x for all i = 1, ..., M 


(S^R 


which is equivalent to 


R c (v) 


gf it d x ~ 


R e (v) 


g d 3 x for all v e (0,1) 


where R c (v) is the region above the fl u =v surface (the subscript 
denotes compliment to distinguish R c (v) from R(v) in Section 
4.5). If the fitting function has the form (5.8), steps similar 
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to those that produced (4.26) can be used to write the above 
equation, for the p-type substrate, as 


B'(v) - B'(l) = q (1 + D h /D e ) (R 0 /V t ) 


g d 3 x . 
. R c< v )' 


( 6 . 1 ) 


The above equation is used to define the best fit for the p— type 
substrate with arbitrary geometry and arbitrary g. Interchange D h 
and D e for the n-type substrate. Integrating this equation with 
respect to v solves for B. The integration constant and B' (1) are 
both determined by the two endpoint conditions B(0)=B(1)=0. 

For the special case of the circular disk in Figure 6.9, fi u is 
well known and the fi u =v surface is seen in the figure as an 
ellipse having the equation 


r 2 sin 2 (7rv/2) + z 2 tan 2 (7rv/2) = r D 2 (equation of n u =v surface). 


Omitting the argument 7rv/2 from the trigonometric functions for 
notational brevity, the integral of g can be written as 


g d 3 x = 2n 

J R C ( V > 


*r D esc 

J 0 


'[r D 2 cot 2 

J 0 


r 2 cos 2 ] 1 /2 

g dz r dr 


which integrates in z first. To integrate in r first, it is 
convenient to make the change in variables w=r 2 sin 2 (ttv/2 ) and 
write the integral as 


g d 3 x = (fr/sin 2 ) 

*r D cot f 

. R c ( v ) J 

0 


z 2 tan 2 

g dw dz . 


The above equations apply to arbitrary g. Specializing to the 
case where g=g Q inside the cylinder and g=0 outside, the integral 
becomes 
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( ff 9o r D 2 ) _1 


g d 3 x = L 

Jr c < v > 


if cos 2 > [(L 2 /2r D 2 ) 2 + L 2 /r D 2 ] 1/2 - L 2 /2r D 2 (6.2a) 


( ff g 0 r D ) 


2x-l 


lR r (v) 


if L 2 /(L 2 +r D 2 ) < cos 2 


d 3 x = L/sin 2 - (2/3) r D cos 4 /sin 3 

- (1/3) (L 3 /r D 2 ) / cos 2 

< [(L 2 /2r D 2 ) 2 + L 2 /r D 2 ) 1/2 -L 2 /2r D 2 (6.2b) 


( w r D ) 


2 \ -1 


g d 3 x = (2/3) (r D /sin) [cos + cos J /(l + cos)] 


lR^(v) 


if cos 2 < L 2 / (L 2 + r D 2 ) 


(6.2c) 


Equation (6.2b) applies when v satisfies the condition that the 
n u = v surface intersects the cylinder wall and lower end, i.e. , 
only the lower cylinder "corners" (seen as corners in Figure 6.9) 
are excluded from R c (v) . The corners contain a small amount of 
carrier generation and there is no need to retain such complexity 
for such an unimportant v interval. Therefore (6.2a) will be used 
over the extended interval cos 2 >L 2 / (L 2 +r D 2 ) . This produces a 
slight discontinuity in the integral of g, equivalent to redis- 
tributing the generation so that the generation in the corners is 
placed on a surface. The total generation within a region that 
completely contains the cylinder is not affected by this redis- 
tribution. Substituting this simplified version of (6.2) into 
(6.1), integrating to solve for B, and replacing L with L-W 
produces the final result 


94 


C-2 



L 0 ■ L - W, L 2 s (L 0 2 + r D 2 ) 1 / 2 , L 2 = Lq/I^, L 3 ^ 

Cq = q ( 1 + D h /D e ) (R 0 /V t ) g 0 r D 2 = 2A D R Q (a Q /N A ) (g 0 /D*) 

C 3 = C 0 L 0' C 2 s (2/3) Cq Tjj 

C 3 = arccos(L 2 ) + (3C 2 /7 r) ln(l/L 3 ) 

+ (3C 2 /tt) In (1 + L 2 ) - (C 2 L 2 /tt) [2 + 1/(1 + L 2 )] 

C v = COS ( 77 v/2) , S v = sin(7T v/2) 

If 0 < C v < L 2 then: 

B'(v) = (C 2 /S v ) (C v + C v 3 /(1 + C v ) ] - C 3 

6(v) = (2 C^/tt) arccos(L 2 ) + (3C 2 /7r) ln(S v /L 3 ) 

+ (3C 2 /7t) In [ 1 + (L 2 - C v )/(1 + C v ) ] - C 3 v 

- (C 2 /7T) [2(1 + C v ) + 1/(1 + L 2 )] [(L 2 - C v )/(1 + C v )] 

If L 2 < C v < 1 then: 

B' (v) = C-l - C 3 
B(v) = (C 3 - C 3 ) v 


95 



which applies to the p-type substrate. For the n-type substrate, 
interchange D h with D e and replace N A with N D in the C Q equation. 

The ambipolar diffusion current G 2 is related to B'(l) via 
(4.26b), and the above results relate B' (1) to C 3 which contains 
the term ln(l/L 3 ). This term becomes singular in the limit as I/-« ° 
(because L 3 -»0) , so the diffusion current has a logarithmic singu- 
larity as the cylinder length is increased without bound in this 
idealized geometry. In a real device, the finite diffusion length 
and/or device dimensions will limit the current if L is too 
large . 

The meaning of a "wide” HRR is interesting. The minority carri- 
er density is negligibly small in a region that extends to infin- 
ity in this idealized geometry . If spatial distance defines HRR 
width, there will always be an infinitely wide HRR. But if a wide 
HRR is to be associated with saturation, we must use something 
other than spatial distance to measure HRR width. Saturation 
occurs when nearly all carrier— modulated substrate resistance is 
in the HRR. A wide HRR will be interpreted to mean that the HRR 
contains nearly all substrate resistance. Because of spreading 
effects, a region that extends to infinity need not have much 
resistance. A wide HRR is not easy to recognize from a plot of 
carrier density versus rectangular coordinates. But it might be 
recognized from a plot that shows how much n u drop is across the 
HRR. This is because a significant n u drop implies that the HRR 
contains a significant fraction of the equilibrium resistance. 
But if the HRR contains a sizable fraction of the equilibrium 
resistance, then it will contain nearly all of the carrier— roodu 
lated resistance (because of the comparatively small carrier- 
modulated conductivity inside of the HRR) . It is therefore most 
informative to plot carrier density against the coordinate v, 
which is related to the spatial coordinates by v=n u (x). The 
actual carrier density will not be a function of v alone, but the 
surface average density on the n u =v surface can be plotted 
against v. The model— predicted density given by (5.11) or (5.24) 
is a function of v alone (constant surfaces are constant 

carrier density surfaces) because the actual g is replaced by a 
fit having the form (5.8). The model-predicted density can be 
regarded as a model-predicted surface average density, and can be 
plotted against v using (5.11) or (5.24) with ft u replaced by v. 

To illustrate saturation and a wide HRR, we consider the spe- 
cific problem of the n + /p diode characterized by 
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r D = 5 nm 


L = 10 /im 

g Q = 6.25 x 10 24 /cm 3 -s 


The value selected for g Q is convenient for making a comparison 
with the first one-dimensional example in Section 6.2. The total 
charge generation rate (including the factor of g) divided by MJ 
area is the same for both cases and is 1000 amps/cm 2 . If the 
device current is normalized by dividing by MJ area (which is 
convenient for making comparisons with Figure 6.1), the normal- 
ized current will be 1000 amps/cm 2 when the device is saturated. 

The model-predicted I-V curve for this three-dimensional exam- 
ple is shown in Figure 6.10. A PISCES prediction is also shown. 
PISCES requires finite geometries and the device simulated by 
PISCES has a cylindrical substrate with a 50-jttm radius and a 50- 
/xm length. The vertical wall is reflective and the lower end is 
the electrode. The version of PISCES used here will not accept a 
g that is uniform inside the cylinder and zero outside. A rough 
approximation of a step function of the radial coordinate is the 
function exp(-r 2 /r D 2 ) , which is the radial dependence used in the 
PISCES simulation. A finite grid spacing results in the PISCES- 
calculated total generation rate being different than the actual 
volume integral of g. An adjusted value was assigned to g Q so 
that the model and PISCES calculate the same total generation 
rate. 

Most of the difference between the two curves in Figure 6.10 is 
due to recombination, which is more noticeable in this extended 
geometry than in the 5— /im one— dimensional geometry represented in 
Figure 6.1. The lifetime for Shockley-Reed-Hall (SRH) recombina- 
tion used by PISCES was 1 microsecond. The difference between the 
two curves in Figure 6.10 is not really very large, but we have 
seen better agreement in Figure 6.1. To improve the agreement and 
verify that other model calculations (e.g., the treatment of a 
three-dimensional geometry) are okay, PISCES was run again for 
the same problem, but with Auger recombination calculations 
turned off and the SRH lifetime changed to 10 milliseconds. This 
virtually eliminates recombination from the PISCES calculations. 
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Figure 6.10: Comparison of I-V curve predictions for the three 
dimensional n + /p diode using r D = 5 nm, L = 10 nm, and g Q = 6.2 
x 10 24 cm -3 



The result is shown in Figure 6.11 and the agreement is now 
almost as good as it was in Figure 6.1. It is more difficult to 
control numerical error associated with grid line spacing in 
three dimensions than in one, and a comparison between PISCES- 
calculated electron and hole currents at various locations to the 
generation rate indicates that there are sizable numerical errors 
in the PISCES curve at small V. Furthermore, the ad hoc fit given 
to PISCES for the radial dependence of g will not produce the 
correct G 2 (unlike the fit used by the model which does produce 
the correct G 2 ) • Both errors are in a direction such that the 
agreement in Figure 6.11 would be further improved at small V if 
these errors were eliminated. This indicates that the model is 
okay except for neglecting recombination. 

It is interesting that the model curves in Figures 6.1 and 6.11 
are indistinguishable. This is not an accident. The two geome- 
tries are almost equivalent as far as the boundary value problems 
(when recombination is neglected) are concerned. The 5-nm one- 
dimensional (ID) problem represented by Figure 6.1 can be given 
the same area 7 rr D 2 as the three-dimensional (3D) problem by in- 
serting a reflective vertical cylindrical wall. This results in 
L=r D for the ID problem, compared to L=2r D for the 3D problem. 
For the ID problem with uniform g, we have G-^=G 2 . Calculating the 
G's for the 3D problem from B'(0) and 6' (1) , which are calculated 
from the equations listed earlier, we find that (when W~0) 
G 2 «l . 06G^ . The total generation rate is the same for the two 
problems, so both G's are nearly the same for the two problems. 
For the ID problem, we have l/R 0 =7ra 0 r D , which is roughly the same 
as the value 4a Q r D applicable to the 3D problem. The two geome- 
tries are not exactly equivalent, but they are almost equivalent. 
The I-V curves saturate early (at small V) , so they are primarily 
controlled by total generation rate and are insensitive to other 
factors such as geometric effects. The small difference in geome- 
try is not observable in these curves. 

The I-V point at V=1 volt in Figure 6.11 is characterized by 
the following model-predicted parameters: 


DR width (W) = 0.403 pm 
substrate voltage drop (V 2 ) = 1.607 volts 
electron density at DRB (P 2 ) = 5.038 x 10^/cm 
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Figure 6.11: Same as Fig. 6. 10 except that recombination was 

eliminated from the PISCES calculations. 


0 



A = 7.256 x 10 12 /cm 3 . 


These data were used to plot surface average electron density 
versus v, as discussed earlier, and the result is shown in Figure 
6.12. This curve can be compared to the curve in Figure 6.2 by 
noting that v is a linear function of distance from MJ in the 
latter figure, satisfying v=l at the DRB and v=0 at the elec- 
trode. The two geometries are not exactly equivalent and the two 
curves are not identical. But plotting density against v instead 
of spatial distance makes a similarity visible. The two curves 
show about the same HRR width if "width" is measured by v instead 
of spatial distance. 


6.5 Conclusions 


Strong funneling is loosely defined by the condition that the 
DR has collapsed and there is a large substrate voltage drop. 
Strong funneling, saturation, and a wide HRR occur together under 
steady-state conditions. A wide HRR in a three-dimensional geome- 
try may be easiest to recognize if surface average minority 
carrier density on the n u =v surface is plotted against v. When 
strong funneling occurs, most substrate voltage is across the HRR 
which limits the current. Because of the high resistance in a 
wide HRR, the current under saturation conditions need not be 
much larger than under nonsaturation conditions. If some parame- 
ter (e.g., the photogeneration rate or a device dimension) is 
varied, the transition between nonsaturation and saturation can 
appear very gradual when looking at terminal currents. Although 
still continuous, the transition appears more abrupt when looking 
at substrate voltage drop or HRR width, and onset conditions can 
be reasonably well defined. The presence of a wide HRR implies 
that the ambipolar diffusion equation fails to provide a good 
approximation for the carrier density function. This equation 
might be used in the AR if boundary conditions are modified to 
account for the presence of the ARB, but a better approximation 
was provided for quantitative estimates. Another observation is 
that carriers need not be generated inside of a DR to collapse 
the DR. Strong funneling can be induced at a distance, i.e., by 
carriers generated outside of the DR. 
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Figure 6.12: Model-predicted surface average electron density on 
the n u = v surface, plotted against v for the three-dimensional 
n + /p diode. 
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A necessary condition for saturation (or DR collapse) was 
derived in terms of ambipolar diffusion currents and is a state- 
ment regarding the spatial distribution of photogeneration. The 
necessary condition is satisfied if all generation is confined to 
a region sufficiently close to the MJ. The condition is not 
satisfied, and a wide HRR cannot form, if the generation is too 
strong at locations too close to the electrode. A saturated 
condition can be, changed to a nonsaturated condition by adding 
additional generation if the addition is near the electrode. This 
does not imply that additional generation can decrease the 
reverse current. The current under nonsaturation conditions can 
exceed the current under saturation conditions if the former case 
is produced by a larger generation rate. If the necessary condi- 
tion for saturation is not satisfied, the DR will not collapse 
even if the carrier density greatly exceeds the doping density in 
the substrate and at the DRB. The condition is more difficult to 
satisfy for the p + /n diode than for the n + /p diode. Compared to 
the n + /p diode, saturation of the p + /n diode requires that gener- 
ation be closer to the MJ. A spatially uniform generation rate 
will not saturate a one-dimensional p + /n diode unless the DR 
width is sufficiently large compared to the substrate thickness. 
When both diode types are operated under similar conditions and 
neither saturates, either can have the larger current. 

The motivation for this steady-state analysis is to obtain 
physical and mathematical guidance for a future transient analy- 
sis. It may be appropriate to point out some of the similarities 
between the steady— state and transient problems. The discussion 
below refers to the transient problem in which funneling is 
induced by an ion that produces a track of free carriers in the 
DR and/or substrate. 

Strong funneling can still be defined by the condition that the 
DR has collapsed and there is a large substrate voltage drop. The 
transient analog of saturation is that the minority carrier 
current at the electrode is negligible at the time of interest. 
It is reasonable to expect this condition to accompany strong 
funneling, and this has been seen in PISCES transient simulation 
results . 

Steady-state funneling can be induced at a distance and there 
is little distinction between carriers generated within the DR 
and carriers generated outside of but close to the DR. Transient 
simulation results show that transient funneling can also be 
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induced at a distance, i.e., not requiring a direct DR hit. 
Carriers must first diffuse to the DR to get the funneling proc- 
ess started. Once there, the DR partially collapses and a sub- 
strate electric field is created. This field drives more minori- 
ty carriers to the DR and the funneling process becomes self- 
sustaining, until the track is sufficiently diminished for the DR 
to recover. Furthermore, carriers generated within the DR do not 
have a special significance. We might expect these carriers to be 
collected much faster than those outside of the DR, because of 
the strong DR electric field, so that these carriers are distin- 
guishable from the others in terms of charge collection time. But 
this is not the case. After the ion hit, carriers initially in 
the DR are separated and driven out as the DR simultaneously 
completely or partially collapses, but charge collection at the 
device terminals does not respond fast enough to be significantly 
affected by this initial charge separation. Terminal currents 
that contribute most to collected charge are seen while the DR is 
recovering, with these carriers now outside of the DR and re- 
sponding to drift and diffusion just like all other nearby carri- 
ers. A plot of collected charge (which is the time integral of 
terminal current and hides current "blips" that negligibly con- 
tribute to collected charge) versus time is very smooth and has 
no demarcation that distinguishes one group of carriers from 
another. "Prompt charge," discussed in some of the older litera- 
ture on single event effects, is not well defined. 

If the track is not long enough to reach the electrode and 
strong funneling. is occurring, it is reasonable to expect an HRR 
below the track, with an electric field inhibiting the downward 
flow of minority carriers. Such an HRR should be depleted of 
minority carriers and (because of quasi-neutrality) depleted of 
excess majority carriers. The visual impression is an absence of 
downward diffusion of the track, and this is seen in transient 
simulation results. 

An interesting case occurs when the track is long enough to 
reach the electrode. Transient simulations were run for two track 
lengths in a reverse-biased n + /p diode. In both cases the DR is 
circular with a 5-/xm radius and was 100 /im above the electrode 
plane. There was a reflective vertical cylindrical wall with a 
50-/im radius. The ion tracks were perpendicular to the device and 
centered on the DR. Both tracks had the same density but one was 
35 /urn long while the other reached the electrode. The collected 
charge versus time curves were almost the same for the two cases. 
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At very early times (< 0.2 ns) after the ion hit, the current 
produced by the long track was a factor of two to three larger 
than that for the short track. At later times (> 0.6 ns) , the 
currents were nearly the same. The time over which the currents 
significantly differed was so short that little charge was col- 
lected during this time, so the collected charge versus time 
curves were nearly the same. At 0.6 ns, the currents were nearly 
the same even though neither track was significantly diminished. 
This situation is similar to the steady-state situation in which 
greatly different HRR widths accompany greatly different sub- 
strate voltage drops to produce nearly the same currents. For the 
long-track transient case, an HRR cannot form until the lower end 
of the track has been cleared away. At 0.6 ns, there has not yet 
been time for much of the track to be cleared away, and the HRR 
is narrow. But because of the small substrate resistance, the DR 
resists collapsing. The DR was only mildly collapsed for the 
long— track case, compared to a greatly collapsed DR for the 
short— track case. The narrow HRR seen in the long-track case was 
accompanied by a correspondingly small substrate voltage drop, 
which produces nearly the same current as the short— track case. 
It is interesting that the long track is less able to maintain a 
collapsed DR than the short track. The spatial distribution 
condition, necessary for steady-state saturation, may have a 
transient analog. 
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APPENDIX A: THE DR EQUATIONS 


A1 Introduction 


The DR equations derived in Reference [2] are very messy. They 
are summarized here for the n"*"/p and p"*"/n junctions, using mostly 
Reference [2] notation, and then simplified. The following nota- 
tion is used: 

U u n = potentials at p-side and n-side DRBs, respectively, 
relative to potential at metallurgical junction 
n p, pp = carrier densities at p-side DRB 

n n , p n = carrier densities at n-side DRB 

V DR = Potential at n-side DRB relative to potential at p-side 
DRB 

H n d = p-side and n-side doping densities, respectively 

j^, j e = scalar current densities at the DRB on the lightly 

doped side; these scalars are positive when currents 
are directed from the n-side towards the p-side 
W = DR width 

S Q = unit step function (5 o (x)=0 if x<0 and 5 0 (x)=l if x>0) 

(i o h , n Q e = low field mobilities 
v = saturation velocity 
V T = thermal voltage 
q = elementary charge 
e = dielectric constant 

a e , a h = 1/ (q V T M 0 , e > and 1 t ^ V T ^o,h> ' respectively 
V T b = 1/ (q v) 

A d = surface area of DRB on lightly doped side 
g = photogeneration rate function 


A2 The n+/p Junction 

The n + /p DR equations were originally listed in Reference [2] 


as 


|Up| = V DR - V T [1 - e VDR/VTj 


(Al) 


[N A +V T b 5 0 (j e )+a e 8 0 ( j e ) V T W/|U p |]W 2 = 2e|U p |/q (A2) 
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n p 2 +[(1/2 )N a -N d e ~ VDR / VT -(3/2)V T b * 0 (j e ) ]n p 

-(1/2) [N A -V T b 5 D (j e )][V T b 5 o (j e )+N 0 e " VDR / VT ) 

-(1/4) [a e 5 0 (j e )] 2 (V T /|U p |)W 2 = 0 . (A3) 


Another equation is needed for which depends on the type of 
conditions considered (steady-state or transient, with or without 
photogeneration in the DR) . The simplest case is steady-state 
with no photogeneration in the DR. For this case, it is often an 
adequate approximation to use j^«0 (a higher approximation is 
available [6] for use when this simple approximation is not 
adequate) . To treat the steady-state case with photogeneration in 
the DR, we use the approximation that the hole current is negli- 
gible in the n + region adjacent to the metallurgical junction 
(MJ) . Then the rate that holes move out of the DR through the p- 
side DRB is simply the photogeneration rate within the DR. Paying 
attention to the sign convention, the equation for j h is 


j h = (q/A D ) 


g d 3 x 

. DR 


q W g D 


(A4) 


where g is approximated by a constant value g D within the DR. The 
complete list of n + /p DR equations consists of (Al) through (A4) . 

When used to solve for.W in terms of the other parameters, (A2) 
is a cubic equation. An exact analytic solution is available, but 
messy, and a simple approximation is more useful. To derive an 
approximation, we first simplify the notation in (A2) by defining 


U - |U p | 
s S W/U 1 / 2 

S Q = (26/q) 1 / 2 [N a + V T b j e ]" 1/2 
K = (q/2e) V T a e j e . 
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We temporarily assume j e >0 so that £ 0 (j e )=l. The results applica- 
ble to j e <0 are trivial to derive and will be listed later. In 
this new notation, (A2) can be written as 


[1/S D 2 + K S/U 1 / 2 ] S 2 = 1 


(A5) 


with S replacing W as the quantity to be solved. It is evident 
from ( A5 ) that for fixed S Q and K, S has the asymptotic forms 
given by 


S -» U 1 / 6 /K 1 / 3 as U -» 0 (A6a) 

S S Q as U -+ oo . (A6b) 

Write (A5) as 

S = [1 - K S 3 /U 1 / 2 ] 1 / 2 S D . (A7 ) 


The strategy is to replace the radical in (A7) with some approxi- 
mation that makes the equation easy to solve for S. The asymp- 
totic forms (A6) show that KS 3 /!! 1 / 2 varies between 0 and 1, so we 
look for an approximation for [1-x 3 ] 1 / 2 that is accurate for 
0<x<l. A particular approximation is 

[1 - x 3 ] 1 / 2 « [1 - x v6 ] 1 / V6 


so that (A7) becomes 


S/S Q * [1 - K V6 / 3 u-vve s V6 3 1/V6 


which can be solved for S with the result 
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s * [S 0 " V6 + K V6 / 3 u-i/vej-i/ve _ 


Going back to the original notation, the equation is written as 

W * (26/q) 1 / 2 lUpl 1 / 2 [(N a + V T b j e ) V6/2 

+ (2e/q) 1 / v6 (V T a e j e ) V6/3 |U p |" 1/V6 r 1/V6 • (A8) 


Numerical calculations will show that (A8) is a very accurate 
approximation of (A2) . 

The original equation (A2) and the approximation (A8) both 
predict the same large | Up | behavior of W, which is 

W « (2e/q) 1/2 | U p I 1 / 2 [N a + V T b j e ]“ 1/2 for large |U p |. (A9) 

It is interesting that (A9) can also be derived by assuming 
velocity saturation in the DR (see any derivation of the Kirk 
effect) . It should be expected that this assumption will lead to 
(A9) because velocity saturation is accompanied by large |U p |. 
But the approximation breaks down for smaller | Up | and it is 
necessary to use the more accurate approximation (A8) . 

Equation (A3) can be greatly simplified, with a small accuracy 
penalty, 1 by replacing W with the large |U p | form given in (A9) . 
The resulting equation can be solved for n p in terms of Vp R (and 
j e ) but the equation is easier to write if V DR is solved in terms 
of n p instead of vice-versa. Solving for the exponential function 
gives 


1 . it is not clear that there is always an accuracy penalty, 
because (A3) is also only an approximation. Table 4.1 of Refer- 
ence [2] compares n_ calculated from (A3) to the values calculat- 

ed from a computer ^simulation for a few special cases. It turns 
out that, for these special cases, (A10) agrees better with the 

computer predictions than (A3) . 
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exp(- V DR /V T ) ~ Nq 1 [Hp - V T b j e ] 

- N d -1 (V T e/q) (a e j e ) 2 [N A 

+ V T b j e ] _1 [N a - V T b j e + 2n p ] 


-1 


(A10) 


This equation restricts the possible values of n p because the 
left side cannot be negative. The allowed values are bounded 
below by the asymptotic (large V DR ) limit, which is the largest 
value that makes the right side of (A10) zero. 

The next approximation is to use | U p | *V DR in (A8). In the 
original derivation, under steady-state conditions with no carri- 
er generation in the DR, there was no distinction between e 
electron current at the DRB and at the MJ. But there is a dis- 
tinction when carriers are generated in the DR, and it is best to 
use the electron current at the MJ in the DR equations Neglect- 
ing the hole current at the MJ, the electron current at the MJ 
the total current at the MJ, which is the total current at the 
DRB. The final modification to the DR equations is to replace j 

with j T ^j h +j e * The final results for the n /P DR ' including th ° Se 
applicable to j T <0, are 


W = (26/q) 1 / 2 Vdr 1 ^ 2 t( N A + V T b 

+ (2e/q) 1 /V6 ,v T a e J T )«/ 3 V6 , -1/V6 if j T > 0 


w = [ (2e/q) Vdr/N ^^] 1 / 2 if j T < 0 


exp (- V DR /V T ) = N d _1 [n p - V T b j T ] 

- N d _1 (V t e/q) (a e j T ) 2 [N A 

+ Vrpb jrp] -1 [N a - V T b j T + 2n p ] 1 if j T > 0 


exp (- V DR /V T ) = N d _1 n p if j T < 0 
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= (q/A D ) 


g d 3 x . 

Jdr 


A3 The p + /n Junction 

The analogous equations for the p + /n DR are 

W = (26/q) 1 / 2 V^ 1 / 2 [(N d + V T b j T ) V6 / 2 

+ (2e/q) 1 / V6 (V T a h j T ) v6 / 3 V^’ 1 ^ 6 ] - 1/V6 if ^ 

W = [ (2e/q) V^/Nd] 1 / 2 if j T < 0 

exp(- V DR /V T ) = H*" 1 [p n - V T b j T ] 

- N a _1 (V t e/q) (a h j T ) 2 [tf D 

+ V T b j^" 1 [N d - V T b j T + 2p n ]’ ] 

exp(- V DR /V T ) = N a 1 p n if j T < 0 


j e = (q/A D ) 


g d 3 x 
DR 


= q W g D . 


> 0 


if j T > 0 
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APPENDIX B: THE SPECIAL FUNCTION H 


B1 Introduction 

Properties of the function H are discussed and a FORTRAN sub- 
program is provided for numerical evaluation. The subprogram can 
be appended to any FORTRAN driver code, allowing the code to call 
the function H as it would call any built-in function. Readers 
not interested in the analytical properties of H can read this 
introduction and then skip to Section B12 on page 134, where the 
subroutine can be found. 

H is loosely defined by the equation 


H(Z 1 ,Z 2 ) = E if and only if e 1/,E = (E - Z 1 )/(E - Z 2 ) • (Bl) 

It is required that either no argument is negative or no argument 
is positive. It is also required that 


1 + Z^ — Z 2 ^ 0 


(B2 ) 


although it is not obvious from a casual inspection of (Bl) why 
(B2) is necessary. The function H has some subtle properties 
requiring a careful analysis, and even the definition has not yet 
been made sufficiently rigorous. Section B2 shows that (Bl) 
sometimes makes sense, i.e., that there exists a unique value for 
H(Z 1# Z 2 ) satisfying (Bl) if Z 1 and Z 2 are suitably restricted. 
Sections B3 through B7 derive bounds, some of which are used in 
Section B8 to take some limits. These limits define H at some 
points that were excluded in Section B2. Nonnegative arguments 
are assumed until Section BIO, which includes nonpositive argu- 
ments. Asymptotic forms are listed in Section B9 , and a suggested 
algorithm for numerical evaluation is given in Section Bll. A 
FORTRAN subprogram, using the suggested algorithm, is listed in 
Section B12. 
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B2 Definition of H(Z1,Z2) when Z1>0, Z2>0, Zl=fZ2, and 1+Zl-Z2=f0 


It is shown in this section that (Bl) makes sense if 


^ 0 / Z 2 > 0 , ^1 4 s ^2' 1 + Z x - Z 2 + 0 (B3 ) 


i.e., that there is a unique E satisfying 


e 1 / E = (E - Z 1 )/(E - Z 2 ) (B4 ) 

if ( B3 ) is satisfied. Note that if we allowed the exponential 
function to have an infinite argument and if Z 1 =0, we would call 
E=0“ a solution to (B4) , where the superscript means that E is on 
the negative side of zero (or 1/E = -«) . Also, if we allowed E to 
be infinite, we would call E=°° a solution. Such cases are not al- 
lowed and (B4 ) does not make sense if E is zero or infinite. 
Existence of a unique E satisfying (B4) means that there is a 
unique nonzero finite E satisfying (B4) . The existence and 
uniqueness proof consists of two steps. The first step proves the 
existence and uniqueness of and X 2 satisfying the three condi- 
tions 


(1 - X x ) e X1 = (1 - X 2 ) e X2 (B5a) 

Z 1 x 2 - Z 2 x i = 1 + z i “ Z 2 (B5b) 

X ± + X 2 . (B5c) 


The second step uses the existence and uniqueness of the X's to 
prove the existence and uniqueness of E satisfying (B4) . 

Before carrying out the first step, it is necessary to define 
and establish some properties of a particular function g (not to 
be confused with the generation rate function) . We start with the 
function f defined by 
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(B6) 


f(X) = (1 - X) e X . 

Differentiating gives f'(X)=-Xe x , so f is strictly increasing 
when X is negative and strictly decreasing when X is positive^ 
Therefore, f is invertible on each branch, i.e., there is an f^ 
and f 2 -1 satisfying 


f 1 ~ 1 (f(X)) = X 

if X 

< 0 

(B7a) 

f 2 _1 ( f (X) ) = x 

if 0 < 

X < 1 . 

(B7b) 

mapping properties are 
f : (-«,0) 4 * (0,1) , 

f ~ 1 • 

r l 

(0,1) « (—oo, 0) 

(B8a) 

f : (0,1) ~ (0,1) , 

f — ^ * 
r 2 

(0,1) «■ (0,1) 

(B8b) 


where «♦ means that the mapping is one-to-one and onto. The func- 
tion f maps the two intervals (-°o,0) and (0,1) onto the same 
target set (0,1), which is the domain of both inverses. Note that 
both inverses are right inverses, i.e.. 


f(f 1 ~ 1 (Y)) = Y and f(f 2 1 ( Y )) =Y if 0 < Y < 1 . (B9) 


The function g is defined by 

g ( 0 ) s 0 (B10a) 

g (X) s f 2 -1 (f (X) ) if X < 0 (BlOb) 

g (X) s f 1 " 1 (f(X)) if 0 < X < 1 (BlOc) 
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and (B9) together with f(0)=l gives 


f (g (X) ) = f(X) if X < 1 . (Bll) 

Using (B8) and (BIO) , we find that g has the mapping properties 

g : (-°°/0) (0,1) (B12a) 

g : (0,1) ** (-«>,0) . (B12b) 


The function g is easiest to visualize from (Bll) and (B12) . 
Given that X<1 and X={=0, we can think of g(X) as "the other X 
producing the same f(X)." In other words, f(g(X))=f(X) but 
g(X)4=X. In fact. 


if X < 1 and X 4= 0/ then g(X) 4= 0 and g( x ) is negative 

(positive) if and only if X is positive (negative) . (B13) 


By combining (BlOa) with (B12) , we get 


g : (~°°, 1) ~ (-°°,i) 


(B14) 


Note that f is strictly increasing on (-°°,0) and decreasing on 
(0,1), so f,” 1 is increasing and f 2 -1 is decreasing. From (B10), 
we conclude that g is decreasing on (-°°,0) and on (0,1). But g is 
continuous at X=0 and we conclude 


g is strictly decreasing on (-»,i) 


(B15) 


Having established some properties of g, we can now show that 
there exists a unique X-^ and X 2 satisfying (B5) . Note that (B5a) 
and (B5c) together imply that X^ and X 2 are both less than 1. 
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This can be shown by contradiction. Assume X 2 >1. Then the right 
side of (B5a) is negative or zero. No negative can make the 
left side negative or zero, so X 2 must be positive or zero. But 
invertibility of f on [0,») contradicts (B5c) . Therefore, (B5) 
implies that X 2 <1. Similarly, X 1 <1. Using the definition of g, we 
can write (B5a) as X 2 =g(X 1 ) . Using this to eliminate X 2 in (B5b) , 
we find that (B5) implies that 

Zj. g(X x ) - Z 2 X 1 = 1 + 7, 1 - Z 2 (B16a) 
x 2 = g^) . (B16b) 


Conversely, (B16) implies (B5) . The equation (B16b) implies (B5a) 
and X 1 ,X 2 <1, and the two equations in (B16) imply (B5b) . Further- 
more, X^=|=0 because the right side of (B16a) is not zero (by 
assumption (B3)). Using (B13) , we conclude (B5c) . Therefore, (B5) 
and (B16 ) are equivalent. To show that there is a unique X 1 and 
X 2 satisfying (B5) , it suffices to show that there is a unique X x 
and X 2 satisfying (B16) . Note that g is strictly decreasing and 
(by assumption (B3)), Z x >0 and Z 2 >0. Therefore the left side of 
(B16a) , regarded as a function of X 1# is strictly decreasing. 
Therefore an Xj^ satisfying (B16a) is unique if it exists. The 
left side maps (-°°,1) onto (- 00 , 00 ) if Z 1 >0. If Z^O, the left side 
maps (-“,1) onto (-Z 2 ,oo). in either case, the target set includes 
the point 1+Z 1 ~Z 2 . Therefore there is a unique X 1 satisfying 
(B16a) . There is a unique X 2 satisfying (B16b) and this completes 
the proof of existence and uniqueness for the X's satisfying 
(B5) . 

We next prove the existence and uniqueness of E satisfying 
(B4). To prove existence, start with the X^ and X 2 satisfying 
(B5) and let 

E = 1/(X 2 - X ± ) . (B17) 

X^ and X 2 exist (are finite) so E=j=0. Furthermore, X 2 =t : X 2 so E 
given by (B17) exists (is finite). By assumption (B3) , Z 1 =}=Z 2 , so 
(B5b) and (B17) can be used to solve for the X's in terms of E 
with the result 
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[Zi - 

(1 + z 2 

- Z 2 ) E]/[(Z 2 - Z x ) 

E] 

(B18a) 

[Z 2 - 

(1 + z x 

- Z 2 ) E]/[(Z 2 - z x ) 

E] . 

(B18b) 


Substituting (B18) into (B5a) shows that E satisfies (B4) , which 
establishes existence of a solution to (B4) . Uniqueness is proven 
by reversing these steps. Let E satisfy (B4) and define X^ and X 2 
by (B18 ) . Using (B18) together with (B4) shows that the X's 
satisfy (B5) , implying that the X's are unique. But E is also 
related to the X's by (B17) , implying that E is unique. 

This completes the proof of existence and uniqueness of E. We 
define H(Z lf Z 2 ) to be this E, so it is now defined for all Z ^ and 
Z 2 satisfying (B3) . 


B3 Some Inequalities 

The function H has been defined when the Z's satisfy (B3) . Some 
other cases such as Z 2 = Z 2 or Z 2 =0 violate (B3) , and limits will 
be used to define H(Z 1 ,Z 2 ) for those cases. Some bounds for 
H(Z lf Z 2 ) will help to evaluate these limits. The first step, and 
the objective of this section, is to derive bounds for the X's 
satisfying (B5) . These bounds will then be used in the next three 
sections to derive bounds for H(Z-^,Z 2 ). The bounds for the X's 
derived here will also be used by the numerical algorithm in 
Section Bll. It is assumed throughout this section that the Z's 
satisfy (B3) so that E satisfying (B4) is related by (B18) to the 
X's satisfying (B5) . 

It was concluded in the previous section that the X ' s are both 
less than 1 and that X-^ is not zero. By combining (B13) with 
(B16b) , we conclude that X 2 is not zero and the two X's have 
opposite signs. Therefore, one of the X's is negative and the 
other is positive and in the interval (0,1). Equation (B5b) can 
be used to identify which of the two X's is negative, and the 
first pair of inequalities is 


X x < 0 and 0 < X 2 < 1 if 1 + 7, 1 - Z 2 > 0 (B19a) 
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0 < x x < 1 


and 


if 


1 + Z± - Z 2 < 0 . (B19b) 


x 2 < 0 


Other bounds can be 
(B5b) . For example, if 
on X 2 via (B5b) , i.e., 


obtained by substituting (B19) back into 
1+Z 1 -Z 2 >0, then X x <0 which implies a bound 


Zl X 2 < 1 + z x - Z 2 if 1 + z x - Z 2 > 0 . 

Similarly, 

Zl x 2 > 1 + Z x - Z 2 if 1 + Z-L - z 2 < 0 • 


(B20a) 


(B20b) 


To obtain (sometimes) tighter bounds, we need a tool derived by 
differentiating e X -(l + X) to conclude that the **1™”*™” 
minimum at X=0. The expression is larger at any X+0 than at X-0, 


or 


,X 


- (l + X) > o 


if 


X + 0 


(B21) 


Now differentiate the expression e -( 1 +X+X /2) and usej > 
conclude that the expression is strictly increasing. The expres 
sion is larger at any x>0 than at X=o, and smaller at any X<0 
than at X=0. This gives 


e x > 1 + X + X 2 / 2 if X > 0 

e x < 1 + X + X 2 / 2 if X < 0 . 


Now write (B5a) as 


(1 - X x ) e _X2 


-XI 


(B22a) 


(B22b) 


= (1 - X 2 ) e 


(B23 ) 



First assume that X 1 <0, implying that X 2 >0. Applying (B22a) to 
the right side of (B23) , and (B22b) to the left side and re- 
arranging terms gives 


(1 - X x ) (1 - X 2 ) < 1 


(B24 ) 


When deriving (B24) , it was assumed that X 2 is negative and X 2 
positive. Interchanging indices for the other case produces the 
same result, so (B24) applies to all cases. Note that (B24) can 
be manipulated into 


X 1 / ( 1 “ * 2 ) > - X 2 (equivalent to (B24) ) . (B25) 

To derive another inequality, note that (B5a) defines X 2 as a 
function of X 2 . Differentiating (B5a) and then using (B5a) again 
to eliminate the exponential function gives 


dX 2 /dX 2 - [U - x 2> X 1 3 / [ (1 - X x ) x 2 ] . 


First assume that X 2 < 0 , implying that 0<X 2 <1. Combining (B25) 
with the above equation gives 


dX 2 /dX2 > x 2 “ 1 > “ 1 • 


The direction of the inequality is preserved upon integration if 
the upper integration limit is larger than the lower. Integrating 
from an arbitrary X 2<0 to X 2=0 while using X 2 ^0 as X 2^0 gives 


X2 t x 2 < 0 . 


(B26) 


When deriving (B26) , it was assumed that XjCO. Interchanging 
indices produces the same result, so (B26) applies to all cases. 

that (B26) states that the negative X has the larger abso- 
lute value. 
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The inequalities (B24) and (B26) were derived from (B5a) alone. 
Including (B5b) allows (B26) to be written as 

X 2 < (1 + Z ± - Z 2 )/(Z 1 + z 2 ) 

and (B24) can be written as 

Z x (1 - X 2 ) 2 + (1 - X 2 ) < Z 2 . 

Temporarily assuming that Z-^Q, this inequality can be manipulat- 
ed into 

[(i - x 2 ) + i/^z-l )] 2 < z 2 /z 1 + 1/iZZ ^ 2 . 

The expression in brackets is positive (because X 2 <1) , so taking 
the square root and rearranging terms gives 

X 2 > 2(1 + Z x - Z 2 ) [1 + 2Z 1 + (1 + 4Z X Z 2 ) X ^ 2 ] ~ 1 

which is also valid when Z 1 =0. 

The important results when (B3) applies, excluding inequalities 
that are always implied by others, are listed below for X 2 
(corresponding bounds for X-^ are implied by (B5b) ) : 


X 2 > 0 if 1 + Z ± - Z 2 > 0 (B27a) 

X 2 < 1 (B27b) 

Z x X 2 > 1 + Z x - Z 2 if 1 + Z ± - Z 2 < 0 (B27c) 

X 2 < (1 + Z ± - Z 2 )/(Z 1 + Z 2 ) (B27d) 
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(B27e) 


X 2 > 2(1 + Z x - Z 2 ) [1 + 2Z X + (1 + 4Z X Z 2 ) x / 2 ] -1 


B4 Bounds for Case 1: 0 < Z2 < Z1 

We assume that (B3) is satisfied and derive some bounds for 
H ( Z^ , Z 2 ) , which is E satisfying (B4) . It is convenient to consid- 
er several cases separately. We start with Case 1, defined by the 
condition 


0 < Z 2 < Z^ (defines Case 1) . (B28) 


E can be expressed in terms of X 2 via (B18b) , with the result 

E = [(Z 2 - Z^ X 2 + 1 + Z x - Z 2 ] _1 Z 2 (B29) 


so that bounds for X 2 imply corresponding bounds for E. Using the 
applicable inequalities in (B27) , and paying attention to the 
fact that Z 1 -Z 2 and 1+Z 1 -Z 2 are positive (for Case 1) when re- 
arranging terms, gives 


E < Z 2 

E < (1/2) (Z-l + Z 2 )/(l + Z x - 
E > Z 2 /(l + Z x - Z 2 ) > 0 . 
Another bound can be obtained by writing (B4) 

E = Z 2 - (Z ± - Z 2 )/(e 1 / E -1) 


z 2 ) 


as 


(B30) 

(B31) 
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The expression on the right, regarded as a function of E is 

. , . /'r P v nl so maps upper bounds into lower 

strictly decreasing (if E>0) , so it maps upp . th 

bounds and positive lower bounds into upper bounds. Using th 
upper bound Z 2 , we obtain the lower bound 


E > Z 2 - (Z-l - Z 2 )/(e 1 / Z2 - 1) • 


Still more bounds can be obtained by writing (B4) as 


E = {ln[(E - Z 1 )/(E - Z 2 )]} 


-1 


(B32) 


Because Z,>Z,, the right side of (B32) , regarded as a function of 
E is strictly decreasing on the interval (0,Z 2 ). Therefore th 
right side of (B32) naps lower bounds for E into upper bounds and 
upper bounds into lower bounds, if E ahd the original bounds are 
inthe reguired interval (0,Z 2 ). But E and the lower bound in 
}b 3 ^ are in the reguired interval, and we obtain the new upper 

bound 


E < { In [ (Z-l + 1) /Z 2 ] } 


-1 


The upper bound in (B30) will be in the required interval if 
Z 2 >l/2, and we obtain the new lower bound 


E > {ln[(Z^ + l/2)/(Z 2 ~ 1/2)]} if Z 2 > ^ 

The bounds for H(Z lf Z 2 ) (=E) are summarized below. 


If 0 < z 2 < z x (Case 1) , then: 


h(z 1; z 2 ) < z 2 


(B33a) 


H(Z 1 ,Z 2 ) < (1/2) (Z x + Z 2 )/(l + Z x - Z 2 ) 


(B33b) 
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H(Z 1# Z 2 ) < {ln[ (Z x + l)/z 2 ]} _1 


(B33c) 


H(Zi,Z 2 ) > Z 2 /(l + Z 1 - Z 2 ) > 0 (B33d) 

H(Z lf Z 2 ) > Z 2 - (Z x - Z 2 )/(e 1 / 22 - 1) (B33e) 

H ( Z 1' Z 2> > ^ ln [( z i + 1 / 2 )/( z 2 - 1/2)]} -1 if Z 2 > 1/2 . (B33f ) 

B5 Bounds for Case 2: 0 < Zl < Z2 < Z1 ± l 
Case 2 is defined by the condition 

0 < Z^ < Z 2 < + 1 (defines Case 2) (B34) 

Some of the applicable inequalities in (B27) combined with (B29) 
give 

0 < Z 2 < E 

(1/2) (Z x + Z 2 )/(l + Z 1 - Z 2 ) < E < Z 2 /(l + Z x - Z 2 ) . (B35) 

Another bound is obtained by using (B18) to write (B24) as 

(E - Z ± ) (E - Z 2 ) < (Z 2 - Z ± ) 2 E 2 
or 

(1+Z 1 -Z 2 ) (l-Z^Zj) E 2 < (Z 1 +Z 2 ) E - Z x Z 2 < (Z 1 +Z 2 ) E . 
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Dividing by the positive quantity (1+Z 1 ~Z 2 ) (l-Z-^+Z^ E gives 
E < (Z 1 + Z 2 )/[(l + Z x - Z 2 ) (1 - Z x + Z 2 )] . 

Because Z 2 >Z^, the right side of (B32) , regarded as a function 
of E, is strictly increasing on the interval (Z 2 ,-h»), which con- 
tains E and both bounds in (B35) if Z 2 >l/2. Therefore, the right 
side of (B32) maps the upper bound into an upper bound and (if 
Z 2 >l/2) the lower bound into a lower bound. The new bounds are 
the same as obtained for Case 1. The bounds for H(Z 1 ,Z 2 ) (=E) are 
summarized below. 


If 0 < Z 1 < Z 2 < Zj + 1 (Case 2), then: 


H(Z 1 ,Z 2 ) < Z 2 /(l + Z 1 - Z 2 ) (B36a) 

H(Z 1 ,Z 2 ) < (Z 1 + Z 2 )/[(l + Z 1 - Z 2 ) (1 - Z ± + Z 2 )] (B36b) 

H(Z 1 ,Z 2 ) < {IntfZ! + 1)/Z 2 ]} _1 (B36c) 

H(Z 1# Z 2 ) > z 2 > 0 (B36d) 

H(Z lf Z 2 ) > (1/2) (Z 1 + Z 2 )/(l + Z x - Z 2 ) (B36e) 


H(Z 1 ,Z 2 ) > {lnC^ + l/2)/(Z 2 - 1/2) ] } _1 if Z 2 > 1/2 . (B36f) 


B6 Bounds for Case 3: 1 < Z1 + 1 < Z2 
Case 3 is defined by the condition 

1 < Z-^ + 1 < Z 2 (defines Case 3) (B37) 
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Some of the applicable inequalities in (B27) combined with (B29) 
give 


(1/2) (Z x + Z 2 )/(l + Z 1 - Z 2 ) < E < Z 1 /(l 
Steps similar to those that produced (B36b) 


+ Z ± - Z 2 ) < 0 . (B38 ) 
give 


E < (Z^ + Z 2 ) /[ (1 + z i “ z 2^ (i ~ z l + z 2^ * 

Because Z 2 >Z 1 , the right side of (B32) , regarded as a function 
of E, is strictly increasing on (-a>,0), which contains E and both 
bounds in (B38) if Z 1 >0. Therefore the right side of (B32) maps 
the lower bound into a lower bound and (if Z^>0) the upper bound 
into an upper bound. The bounds for H(Z-^,Z 2 ) (=E) are summarized 
below. 


If 1 < z x + 1 < Z 2 (Case 3), then: 

H(Z 1 ,Z 2 ) < z 1 /(l + Z ± - z 2 ) < 0 (B39a) 

H(Z 1 ,Z 2 ) < (Z 1 + Z 2 ) / ( (1 + Z 1 " Z 2 ) (1 - Z x + Z 2 ) ] < 0 (B39b) 
H(Z 1 ,Z 2 ) < {ln[Z 1 /(Z 2 - l)]}" 1 if Z x > o (B39C) 

H(Z 1# Z 2 ) > (1/2) (Z 1 + Z 2 )/(l + Z ± - Z 2 ) (B39d) 

H(Z 1 ,Z 2 ) > {ln[ (Z^ + l/2)/(Z 2 - 1/2) ] } -1 . (B39e) 
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B7 Some Addition al Bounds for X2 

The bounds for E (-H(Z lr 2 2 )» derived in the last three sections 
are adequate for the intended purpose of determining a f 
selected limits and asymptotic forms (next two sections) . B 
u^ss a limit or asymptotic form is found to apply, numerica 
evaluation of H(Z 1( Z 2 ) will work with the X-s and the number of 
required calculations is reduced by tightening the bounds for X 
It is therefore desirable to use all available information to 
bracket X 2 as tightly as possible. Some of the bound. .for E are 
equivalent to (via (B18b) ) or weaker than the X 2 bounds ™ < >; 

some other E bounds, such as <B33c) , were obtained from sn addi 
tional step and can be used to derive new X 2 bounds. Using 
^B18b), (B33C), (B33e) , (B33f) (B360). <B36f), (B39C), and 

(B39e) provides the following additional bounds. 



if 0 < z 2 < z x 




Z 1 +l/2 

Z 2 -l/2 


if Z^>0, Z 2 >l/2 / Z 2 >Z! 

and 1 +Z^— Z 2 ^0 


X 2 



+ 




Z 1 +l/ 2 

z 2 -i / 2 


if 1/2 < Z 2 < 


l+Zi“Z 2 




if 0 < Z 1 < Z 2 < Z 1+ l 



if 1 < Z-l+1 < Z 2 


127 


x > 1+Z 1 Z 2 + Z 2 1 ~ e" 1 / 22 if 0 < Z 2 < Zl and 

Z 1“ Z 2 Z 2“ Z 1 Z 2 “ Z 1 e -1 / 22 Z 2 < [ln(Z 1 /Z 2 ) ] _1 


B8 Definition of HfZ.Z) and HfZ.cn when Z>0 

The quantities H(Z,Z) and H(Z,0) are not yet defined because 
the arguments violate (B3) . These quantities will be defined, 
with Z>0, by taking limits. A limit of a function of several 
variables can be subtle because a given point can be approached 
along a variety of paths, and the limit is well defined only if 
all possible paths produce the same limit. Fortunately, the 
limits needed here are well defined. 

First consider the limit as (Z lf Z 2 ) approaches (Z,Z) for some 
Z>0. If (Z 1 ,Z 2 ) is sufficiently close to (Z,Z), Case 3 is exclud- 
ed and the bounds (B33a) , (B33d) , (B36a) , and (B36d) imply that, 
no matter what path is followed, we have H(Z 1 , Z 2 )->Z. We define 
H(Z,Z) to be this limit, i.e.. 


H ( Z , Z) = Z if Z > 0 . (B40a) 


Now consider the limit as ( Zl ,Z 2 ) approaches (Z,0) for some 
Z>0. We may assume that Z>0, because (B40) applies if Z=0. But if 
Z>0 and ( Z x , Z 2 ) is sufficiently close to (Z,0), only Case 1 can 
apply. The bounds (B33a) and (B33d) imply that, no matter what 
path is followed, we have H( Z;L ,Z 2 )-o. We define H(Z,0) to be this 
limit, i.e., 


H (Z , 0) = 0 if Z > 0 . 


(B40b) 


The condition Z— 0 was allowed in (B40b) because (B40a) and (B40b) 
are equivalent when Z=0. The quantity H(Z-^,Z 2 ) is now defined for 
all nonnegative Z ^ and Z 2 satisfying 1+Z 1 -Z 2 =|=0. 
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B9 Asymptotic Forms 


Asymptotic forms are approximations for H(Z^,Z 2 ) that become 
exact (in the sense that the relative or fractional error goes to 
zero) in the limit as various combinations of the arguments 
become small or large. Such approximations make the behavior of H 
easier to visualize. They also have computational advantages 
(when applicable) because they are simple. With (B3) assumed, 
asymptotic forms are derived below for: small | Z^_— Z 2 1 # small Z 2 , 
small | 1+Z 1 -Z 2 | , and large Z ± and Z 2 * Applicability tests are 
given in terms of a positive quantity « tol , which is a user 
specified relative error that will be tolerated in the estimate 
of H(Z 1 ,Z 2 )* For example, if an error less than one percent is 
good enough, whether too large or too small, then 5 tol =0.01. 

The first asymptotic form applies when e^= | z i~ z 2 I i s small. A 
sufficiently small excludes Case 3. To insure that Case 3 is 
excluded, it is required that £ 1 <1. For Case 1 conditions, 
e 1 =Z 1 -Z 2 and the bounds (B33a) and (B33d) can be written as 


Z 2 /(l + e x ) < H(Z 1 ,Z 2 ) < Z 2 


Similarly, Case 2 gives 


Z 2 < H(Z 1 ,Z 2 ) < Z 2 /(l - e i) . 


In either case, the approximation H(Z 1 ,Z 2 )~Z 2 has a relative 

error less than e 1# if e ± <l. The approximation has a relative 

error less than ^^ol ^ e l < ^ an< ^ e l < ^tol' i* e *» 

H(Z 1 ,Z 2 ) * Z 2 has relative error less than ^tol if ( B3 ) (B41a) 

applies and |Z 1 -Z 2 | < 1 and | Z 1 ~Z 2 | < ^tol* 

The second asymptotic form applies when Z 2 is small. A suffi- 
ciently small Z 2 excludes Case 3 conditions. A small Z 2 under 
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Case 2 conditions implies a small \Z 1 -Z 2 \ and (B41a) applies. It 
is therefore adequate to consider only Case 1 conditions by 
requiring that 0<Z 2 <Z^. "Small Z 2 " under Case 1 conditions is 
interpreted to mean that e 2 is small, where 

e 2 = [(Z x - Z 2 )/Z 2 ] [e 1 / 22 - l]" 1 . 


The Case 1 bounds (B33a) and (B33e) can be written as 


(1 e 2 ) Z 2 ^ H ( Z^ / Z 2 ) ^ Z 2 • 


In addition to 0<Z 2 <Z 1 , we also require that the Z's satisfy the 
condition e 2 <l. T ^en the above bounds imply that the approxima- 
tion H(Z^,Z 2 )«Z 2 has a relative error less than e 2 /(l-e 2 ). The 
approximation has a relative error less than $tol if the require- 
ments 0<Z 2 <Z 1 , e 2 <l, an< * e 2^ <l ^tol are sa tisfied. Using 
the definition of e 2 to express the e 2 requirements in terms of 
the Z's and noting that one of the inequalities is implied by the 
other, we obtain 


H(Z lf Z 2 ) * Z 2 


has relative error less than ^^ol ^ 

(B3 ) applies and Z 2 <Z 3 and 

Z 2 < {ln[Z 1 /Z 2 + (Z 1 -Z 2 )/(Z 2 * tol )]} 


-1 


(B41b) 


The third asymptotic form applies when the Z's come close to 
the forbidden condition 1+Z,-Z 2 =0, i.e., e 3 s|i+Zj-Z 2 | is small. A 
sufficiently small e 3 excludes Case 1 conditions. We insure that 
Case 1 is excluded by requiring that e 3 <l. For Case 2, e 3 =l+Z 1 ~Z 2 
and (B36b) and (B36e) can be written as 

[(1/2) (Z 1 + Z 2 )/(l + Zj - Z 2 ) 3 < H(Z 1 ,Z 2 ) 

< [(1/2) (Z 1 + Z 2 )/(l + Z ± - Z 2 )]/(l - e 3 /2) . 


Similarly, Case 3 gives 
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[ - (1/2) (Z x + Z 2 ) / (1 + Z 1 - Z 2 )]/(l + C 3 / 2 ) < _ H(Z 1 ,Z 2 ) 

< [ - (1/2) (Z x + Z 2 )/(l + Z ± - Z 2 ) ) . 

In either case, the approximation 

H(Z 1 ,Z 2 ) « (1/2) (Z x + Z 2 )/(l + Zi “ Z 2> 

has a relative error less than e 3 /2, if c 3 <l. The approximation 
has a relative error less than 5-tol e 3 <:L and e 3/ 2 < 5 tol' ^* e *' 

H(Z 1 ,Z 2 ) « (1/2) (Z x + Z 2 )/(l + z i " z 2 > 

has relative error less than ^tol (B41c) 

(B3) applies and 
1 1 +Z-L- Z 2 1 < min{ 1 , 2<S tol } . 

The last asymptotic form applies when the Z's are both large in 
the sense that the upper and lower logarithmic bounds ((B33c) and 
(B33f) for Cases 1 and 2, and (B39c) and (B39e) for Case 3) come 
together. An equivalent statement is that e 4 is small, where 

e 4 a 1 - ln[(Z 1 + 1 ) /Z 2 ]/ln[ (Z 3 + l/2)/(Z 2 - 1/2)] 

if Z 2 > 1/2 and 1 + Z x - Z 2 > 0 

e 4 = ln[Z 1 /(Z 2 - l)]/ln[(Z 1 + l/2)/(Z 2 - 1/2)] - 1 

if Z 1 > 0 and 1 + Z 3 - Z 2 < 0 . 

The logarithmic bounds can be written as 
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{l/l n t ( z i + l/2)/(Z 2 - 1/2)]} < H(Z 1 ,Z 2 ) 

< {l/ln[ (Z x + l/2)/(Z 2 - 1/2) ]}/ (l - e 4 ) 

if Z 1 >0, Z 2 > 1/2, Z 2 Z lf and 1 + Z 1 - Z 2 > 0 

{ - l/ln^Z-L + l/2)/(Z 2 - 1/2) ]} / (1 + e 4 ) < - H(Z 1 ,Z 2 ) 

< { - l/ln[(Z 1 + l/2)/(Z 2 - 1/2)]} 

if Z 1 > 0 and 1 + Z-± - Z 2 < 0 . 

In either case, the approximation 

H(Z 1 ,Z 2 ) * {lnC^ + l/2)/(Z 2 - 1/2) ] } _1 

has a relative error less than e 4 . The relative error will be 
less than $ tol if e 4 <S tol . The fact that ln(a+b) <ln(a) +b/a, when 
a>0 and b>0, can be used to show that if Z 1 >0, Z 2 >l/2, and 1+Z^- 
Z 2 =^0, then 

e 4 < (1/22^ [ (l+Z 1 -Z 2 )/(Z 2 -l/2) ] {ln[ (Z.,+1/2) / (Z 2 ~l/2) ] } _1 

so the relative error is less than $ tol if the Z's satisfy the 
stated conditions and the right side of the above inequality is 
less than $tol* This gives 

H(Z 1 ,Z 2 ) « {lnt^ + l/2)/(Z 2 - 1/2) ] } _1 

has relative error less than if (B41d) 

(B3 ) applies and Z 1 >0, Z 2 >l/2, and 

1 /Z 1 < 2 S tol [(Z 2 -l/2)/(l+Z 1 -Z 2 ) ] ln[(Z 1 +l/2)/(Z 2 -l/2)] . 
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BIO Definition of HfZl.Z2) w hen Z1<0, — Z2 < 0 , — and 1+Z1-Z2={=0 

So far, H(Z^,Z 2 ) has been defined for any nonnegative Z's 
satisfying l+Z 1 -Z 2 +0. The definition is either (Bl) or limits 
derived from it. A similar definition can be given when the 
arguments are nonpositive. For this case, we can manipulate (Bl) 

into 


e l/(-E) = t (-E) - (-Z 2 )]/[(-E) - (-Zi)] 


which has the same solutions and limits as the original (Bl) , 
except for a change in symbols; -E replaces E, -Z 2 replaces Z 1# 
and -Z ± replaces Z 2 . The problem case is l+(-Z 2 ) -(-Z 1 )=0, so even 
for nonpositive arguments the problem case is still 1+Z 1 -Z 2 =0. We 
therefore define 


H(Z 1 ,Z 2 ) ■ - H(-Z 2 ,-Z 1 ) if Z,<0, Z 2 <0, 

so H(Z 1 ,Z 2 ) is now defined for all Z's 
1+Z^-Z 2 =^0 . 


and 1+Z^-Z 2 4=0 (B42) 

satisfying Z 1 Z 2 >0 and 


Bll A Numerical Algorithm 

This section suggests one possible algorithm for numerical 
evaluation of H(Z 1 ,Z 2 ). Because of (B42) , it is sufficiently 
general to confine our attention to nonnegative arguments. It is 
assumed below that the Z's are nonnegative and satisfy l+Z 1 -Z 2 +0. 

The first step is to determine the applicability of the special 
cases (B40) and asymptotic forms (B41) in the order listed. Use 
the first case that was found to apply. If none apply, then the 
Z's satisfy (B3) , implying that H(Z 1 ,Z 2 ) is the E satisfying (B4) 
and can be solved via (B29) from the X's satisfying (B5) . The 
equations governing the X ' s can be written as 


X 1 = [Z 1 X 2 - (1 + Z-l - Z 2 )]/Z 2 


(B43a) 
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W = X x - X 2 + ln[ (1 - X^/fl - X 2 )] 


(B43b) 


W = 0 . 


Let X A and X B be, respectively, the largest lower bound and 
smallest upper bound for X 2 that can found in (B27) and in 
Section B7 . The basic idea is to guess at X 2 , then calculate X^ 
from (B43) and W from (B44) . Whether the guess is too large or 
too small depends on the sign of W. The guess is too large if W 
has the same sign as W B , where W B is the value obtained when X 2 
is replaced with X B in (B43) . The guess is too small if W and W B 
have opposite signs, and the guess is correct if W=0. The bisec- 
tion method is used to construct a seguence of lower X 2 bounds 
X A ^> , X A < 2 > , ..., and upper bounds Xg^ 1 )), Xg( 2 )), .... The 
first bounds are X A ( 1 )=X A , Xg^^Xg. For i>l, the i th bounds are 
constructed from the previous bounds by letting be the 
midpoint 


x M (i ~ 1} = (1/2) [X A ( i_1 ) + Xgf 1 " 1 )) . 


Determine the sign of W obtained when X 2 is replaced with this 
midpoint. If W and W B have the same sign, the correct X 2 is 
smaller than the midpoint and the new lower bound is the same as 
the old while the new upper bound is the old midpoint. Similarly, 
if W and W B have opposite signs, the new lower bound is the old 
midpoint and the upper bound is not changed. As the upper and 
lower X 2 bounds come together, the corresponding E bounds (from 
(B29) ) also come together. The bisection is terminated when the E 
bounds are sufficiently close together. 


B12 The Function Subprogram 

The subprogram listed at the end of this discussion can be 
appended to any FORTRAN source code, allowing the code to call 
the function H as it would call any built-in function. The sub- 
program uses the numerical algorithm discussed in the previous 
section. 
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This computer version of H differs from the analytical version 
in that there is a redundant argument Z 3 =Z 1 -Z 2# included to 
improve numerical accuracy. It is desirable for a computer code 
to be able to deal with nearly any case allowed by the mathemati- 
cal theory. One allowed case is that in which Z 3 and Z 2 are 
nearly equal in the sense that the difference Z 1 -Z 2 is a tiny 
fraction of either of the two Z's. If the difference is calculat- 
ed by letting the computer subtract the nearly equal Z's, the 
relative error will be large unless the Z's are passed to H with 
sufficiently high precision and the subtraction performed with 
the same precision. Error is especially disruptive if the differ- 
ence is close to the forbidden value -1, because H is singular 
when 1+Z,-Z 2 =0. An alternative to passing two high-precision 
arguments (and performing high-precision arithmetic) is to pass 
three lower precision arguments with the difference being one of 
the arguments. Of the three arguments, only the two having the 
smallest absolute values really need be passed. But the subpro- 
gram allows Z. and Z 2 to be any pair of numbers in the domain of 
H so it is not known in advance which argument has the largest 
absolute value, and any one of the three can have the smallest 
absolute value. Therefore, all three arguments are passed. 

A tolerance parameter DELTOL is set equal to 10 -4 . This would 
result in H(Z.,Z 2 ,Z 3 ) being calculated with an error of less than 
one part per ten thousand (the intended accuracy) , if machine 
precision was unlimited. An effort was made to manipulate expres- 
sions into forms that do not subtract nearly equal numbers. In 
spite of this effort, machine precision can still limit the 
accuracy in some extreme cases, such as when Z 3 is very close to 
-1 where H is undefined. The intended accuracy is not always 
guaranteed, and will still not be guaranteed even if DELTOL is 
assigned a smaller value. 


C 

C 

C 

C 


This F function"subprogram can be appended to a FORTRAN source 
code allowing the code to call the function H defined in the 
?ex?: This computer version of H differs from the analytical 
c version in that there is a redundant argument Y3-Y1 Y2. A 
C three arguments are passed to insure that the two having the 
c smallest absolute values are represented with the greatest 
C possible numerical precision. Y1 and Y2 cannot have opposite 
C signs, and Y3 cannot be -1. 

C 

DELT0L=1 . OE-4 
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ooooo ooooo o on 


Check for illegal arguments. 

IF (Y1*Y2 .LT. 0. 0) THEN 

WRITE (*,*) 'ERROR: FIRST TWO ARGUMENTS HAVE OPPOSITE SIGNS' 
GO TO 130 
END IF 

IF (Y3.EQ.-1.0) THEN 

WRITE (*,*)' ERROR: THIRD ARGUMENT IS -1» 

GO TO 130 
END IF 

Assign a new value to the argument having the largest absolute 
value if needed to comply with Y3=Y1-Y2, without disturbing 
the other two arguments. 

IF (Y2*Y3 .GE. 0. 0) Y1=Y2+Y3 
IF (Y1*Y3 . LE. 0.0) Y2=Y1-Y3 

If Y1 and/or Y2 are negative, use H(Y1,Y2,Y3)= -H(-Y2,-Y1, Y3) . 
The arguments used will be Zl, Z2, and Z3. Set a flag as a 
reminder to multiply H by -1 if the Z's differ from the Y's. 

IFLAG=0 
Z1=Y1 
Z2=Y2 
Z3=Y3 

IF ( (Y1.LT.0.0) .OR. (Y2.LT.0.0) ) THEN 
Z1=-1.0*Y2 
Z2=-1.0*Y1 
IFLAG=1 
END IF 
C 

C Use special cases or asymptotic forms if applicable. T with 
C or without subscripts is for temporary storage of intermediate 
C results and can represent different quantities in different 
C calculations. 

C 

IF (Z3.EQ.0.0) THEN 
E=Z2 

GO TO 120 
END IF 
C 

IF (Z2.EQ.0.0) THEN 
E=0 . 0 
GO TO 120 
END IF 
C 

AZ3=ABS(Z3) 

IF (AZ3 .LT.DELTOL) THEN 
E=Z2 

GO TO 120 
END IF 
C 

IF (Z3.LE.0.0) GO TO 10 
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T=Z3* (DELTOL+1 . 0) / (Z2*DELTOL) 
Tl=ALOG ( 1 . 0+T ) 

T2=l . O/Tl 

IF (Z2.LT.T2) THEN 
E=Z2 

GO TO 120 
END IF 

10 CONTINUE 
C 

AZD=ABS ( 1 . 0+Z3 ) / 2 . 0 
IF (AZD.LT.DELTOL) THEN 
E=0 . 5* ( Z1+Z2 ) / (1. 0+Z3) 

GO TO 120 
END IF 


IF (Z1.LE.0.0) GO TO 20 
IF (Z2.LT.1.0) GO TO 20 
T= ( Zl+0 . 5 ) / (Z2-0.5) 

Tl=ALOG(T) 

T2=2.0*DELTOL* (Z2-0. 5) *T1/ (1.0+Z3) 
T3=l. 0/T2 

IF (Z1.GT.T3) THEN 
E=1 . O/Tl 
GO TO 120 
END IF 

20 CONTINUE 


C 

C 

C 

C 

C 

C 


If none of the above cases apply, prepare to estimate X2 by 
constructing a lower bound XA and an upper bound XB. Start 
with bounds that always apply and then go through the list to 
see whether the upper bound can be made smaller or the lowe 
bound made larger . 


XB=1. 0 

T=l. 0+4 . 0*Z1*Z2 
T1=SQRT (T) 

T2=1.0/ ( 1 . 0+T1+2 . 0*Z1) 
XA=2 . 0*T2* ( 1 . 0+Z3 ) 

IF (Z3.GT.-1.0) GO TO 30 
IF (Z1.LE.0.0) GO TO 30 
XAN= ( 1 - 0+Z3 ) /Z1 
IF (XAN.GT.XA) XA=XAN 
30 CONTINUE 

XBN= ( 1 . 0+Z3 ) / (Z1+Z2) 

IF (XBN.LT.XB) XB=XBN 

IF (Z3.LT.0.0) GO TO 40 
T= ( 1 • 0+Z3 ) / Z2 
Tl=T-ALOG ( 1 . 0+T) 
XBN=Z2*T1/Z3 
IF (XBN.LT.XB) XB=XBN 
40 CONTINUE 


C 
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IF (Z2.LT.1.0) GO TO 50 
IF (Z3.GT.0.0) GO TO 50 
T= ( 1 . 0+Z3 ) / (Z2-0.5) 

Tl=ALOG ( 1 . 0+T) 

T2=Z2*T1/Z3 
T3= ( 1 . 0+Z3 ) / Z3 
XBN=T3 -T2 

IF (XBN.LT.XB) XB=XBN 
50 CONTINUE 
C 

IF (Z2.LT.1.0) GO TO 60 
IF ( Z 3 . LT .0.0) GO TO 60 
T=(1.0+Z3) / (Z2-0.5) 

Tl=ALOG ( 1 . 0+T) 

T2=Z2*T1/Z3 
T3= ( 1 . 0+Z3 ) / Z3 
XAN=T3-T2 

IF (XAN.GT.XA) XA=XAN 
60 CONTINUE 
C 

IF (Z3.LT.-1.0) GO TO 70 
IF (Z3.GT.0.0) GO TO 70 
T= ( 1 . 0+Z3 ) / Z2 
Tl=T-ALOG(l. 0+T) 

XAN=Z2*T1/Z3 
IF (XAN.GT.XA) XA=XAN 
70 CONTINUE 
C 

IF (Zl.LE.O.O) GO TO 80 
IF (Z3.GT.-1.0) GO TO 80 
IF (Z2.LT.1.5) GO TO 80 
T=(1.0+Z3) / (Z2-1.0) 

Tl=ALOG ( 1 . 0+T) 

T2=Z2*T1/Z3 
T3= ( 1 . 0+Z3 ) / Z3 
XAN=T3-T2 

IF (XAN.GT.XA) XA=XAN 
80 CONTINUE 
C 

IF (Z3.LT.0.0) GO TO 90 
TO=ALOG ( 1 . 0+Z3 / Z2 ) 

T0=1 . 0/T0 

IF (Z2.GE.T0) GO TO 90 
T=EXP (-1.0/Z2) 

Tl= (1 . 0-T) / (Z2-Z1*T) 

T2=Z2*T1/Z3 
T3= ( 1 . 0+Z3 ) / Z3 
XAN=T3-T2 

IF (XAN.GT.XA) XA=XAN 
90 CONTINUE 
C 

C If XA and XB are so close together that numerical error gave 
C XB<=XA, use X2= ( 1/2 ) (XA+XB) and calculate E and skip the 
C bisection loop. 

C 
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noon non 


IF (XB.LE.XA) THEN 
X=0 . 5* (XA+XB) 

T= ( 1 . 0+Z3 ) -Z3*X 
E=Z2/T 
GO TO 120 
END IF 


C 

C Now 


calculate WA and WB and record their signs in SA and SB 


X1=(Z1*XB— (1.0+Z3) ) / Z2 
WB=1 • 0 

IF (XB.LT.1.0) WB=X1- XB+ALOG ( ( 1 . 0— XI ) / ( 1 . 0— XB) ) 

SB=1 . 0 

IF (WB . EQ .0.0) SB=0 . 0 
IF (WB . LT .0.0) SB=-1.0 

Xl= ( Z1*XA- ( 1 . 0+Z3 ) ) /Z2 
WA=X1-XA+AL0G( (1.0-X1) / (1.0-XA) ) 

S A= 1.0 

IF ( WA . EQ .0.0) SA=0 . 0 
IF (WA. LT. 0.0) SA=-1.0 


C 

c 

c 

c 

c 

c 

c 


If XA or XB are so close to the correct solution X2 that 
™ m 2ical error gave XB<=X2 or XA>=X2 SA*SB will be positive 
or zero. If SA*SB= -1, everything is okay and the next block 
of steps can be skipped. Otherwise, determine which of the 
intended bounds is closest to X2 . Set X2 equal to that 
intended bound and calculate E and skip the bisection loop. 


IF ( S A* SB . EQ .-1.0) GO TO 100 
T=(l. 0+Z3) -Z3*XA 

IF (SB*WB. LT. SA*WA) T= ( 1 . 0+Z3 ) -Z3*XB 
E=Z2/T 
GO TO 120 
100 CONTINUE 


NOW start the bisection loop to tighten the X2 bounds. 


110 CONTINUE 

X=0 . 5* (XA+XB) 

Xl= ( Z1*X- ( 1 . 0+Z3 ) ) / Z2 
W=X1-X+AL0G( (1. 0-X1) / (1. 0-X) ) 

S=l. 0 

IF (W.LT.0.0) S=-1.0 
IF ( S*SB . GT .0.0) XB=X 
IF (S*SB . LT. 0.0) XA=X 

If W=0, the solution was found. Calculate E and exit from the 
loop. 


IF (W. EQ. 0.0) THEN 
T= ( 1 • 0+Z3 ) -Z3*X 
E=Z2/T 
GO TO 120 
END IF 
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c 

C EA and EB are the values of E when X2 is replaced with XA and 
C XB respectively. If EA and EB are close enough together, 

C calculate E and exit from the loop. Otherwise, go through the 
C loop again. 

C 

T= ( 1 . 0+Z3 ) -Z3 *XA 
EA=Z2/T 

T= (1 . 0+Z3 ) -Z3*XB 
EB=Z2/T 

E=2 . 0*EA*EB/ (EA+EB) 

T= (EA-EB) / (EA+EB) 

DELTA=ABS (T) 

IF (DELTA. LT. DELTOL) GO TO 120 
GO TO 110 
120 CONTINUE 
H=E 

IF (IFLAG.EQ.l) H=-1.0*E 
130 CONTINUE 
RETURN 
END 
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APPENDIX C: THE SPECIAL FUNCTION F 


The function F is defined by Y=F(X 1 ,X 2 ) if and only if Y satis 
f ies 


Y + (1 - X x ) ln(l + Y/X x ) = X 2 . 


(Cl) 


It is sufficiently general to confine our attention to those 
cases where X 1 is positive and X 2 is positive or zero. X 1 is not 
allowed to be zero. 

The function F is closely related to a particular type of 
inverse of the special function H. If we want to solve (Cl) for 
X n when Y and X 2 are given, the solution can be expressed in 
terms of H. If we want to solve (Cl) for Y with the X’s given, we 
use F. But F is much simpler than H and an approximation for F 
was already listed in the form of the generalized ambipolar 
approximation. The connection is made clear by remembering where 
F first originated. For the p-type substrate with g=0. Section 
3.2 found that 


P + (p 0 /2 - A) ln(l + P/A) = n 


or 


P = (p 0 /2) F (2A/p 0 , 2fl/p 0 ) 


(C2 ) 


The generalized ambipolar approximation is an approximation for 
either side of (C2) . 

Iterations are used to evaluate F, or solve for Y. Iterations 
are performed by manipulating (Cl) into 


Y = f(X 1# X 2 , Y) ( C3 ) 


for some appropriately chosen f, which is not unique. If f is 
chosen well, the sequence of iterates Y ( Y ( > , ... will con- 
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verge to the solution, where Y ( is some initial guess and 
Y (i+l) = f(Xi/ X2> Y (i) } (i = 0t lt # 

Several cases, characterized by the way X 2 and X 2 compare with 
each other, are considered separately. The need for considering 
different cases can be seen from the fact that P given by (C2) 
behaves differently under different conditions. If A>p Q /2 (imply- 
ing that V 2 <0) , the nominal ambipolar approximation is a good 
low-order approximation. If A<p Q /2 (implying that V 2 >0) , there is 
an HRR and an AR, and the behavior of P depends upon which region 
we are examining. We can anticipate that at least three cases 
require separate treatment. It turns out that there are four 
cases, with one corresponding to a transitional region near the 
ARB. 


The different cases will use different f’s in (C3) and differ- 
ent intervals from which the initial guess is to be selected. The 
proof of convergence is fundamentally the same for all cases. The 
basic idea is to find a closed interval, from which the initial 
guess is to be selected, having the property that f, regarded as 
a function of Y, maps this interval into itself. Then show that, 
throughout this interval, the absolute value of the Y derivative 
of f is less than or equal to some number that is strictly less 
than 1 (preferably less than 1/2 so that the iteration will 
converge at least as fast as the bisection method) . The details 
are omitted because they are not difficult. Error estimates are 
also obtained by iteration, but not necessarily the same conver- 
gent iteration that produces progressively better estimates. The 
basic idea is the same for all cases and illustrated for the 
first case considered. 

We start with Case 1 defined by 


X^ > 1 (defines Case 1) . 


This case will be encountered when we want to use (C2) to solve 
for P and there is no HRR (i.e., A>p Q /2 or V 2 <0) . The iteration 
is 
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y(i+l) = x 2 + (x x - 1) In ( 1 + Y^/Xj.) (for Case 1) (C4) 


which converges for any initial guess selected from the interval 
[X 2 ,°°). The suggested initial guess is 

y(°) = X 2 (for Case 1) . 


Convergence of the iteration (C4) can be very slow in theory. In 
practice, Case 1 is accompanied by X-^l and the convergence is 
fast. An error estimate associated with any given iterate is 
obtained by manipulating (Cl) into 


y = Y + [(X 2 + X 1 )/(X 2 + 1)3 [X 2 + (X! - 1) ln(l + Y/X x ) - Y] . 

The actual solution Y and all iterates produced by (C4) are in 
the interval [X 2 ,«). Differentiating shows that the right side of 
the above equation is decreasing in Y on this interval (or con- 
stant if X 1 =l) . Therefore the right side maps iterates that are 
too small into estimates that are too large and vice-versa. The 
correct solution is bracketed by any iterate Y^ 1 ' and its conju- 
gate Yq^ 1 ) defined by 

y c (i) s Y ( i) 

+ [ (X 2 +X-L) / (X 2 +l) ] [X 2 +(X 1 -1) ln(l+Y (l) /X^- Y* 1 *] (for Case 1) . 

The difference between Y ^ ^ and Y c ^^ is a simple error estimate. 
Case 2 is defined by 


0 < X^ < 1 and 

X 2 > 2 + (1 - X x ) In ( 1 + X^-Xj.) (defines Case 2) . 

It can be shown that Case 2 implies that Y>2, or P>P G in (C2) . 
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This case is encountered when there is an AR and HRR (A<p Q /2) and 
we want to use (C2) to solve for P at some point in the AR not 
too close to the ARB. The iteration is 


y( 1+1 ) = + [ (X 2 + X 1 )/(X 2 + 1)] [X 2 

+ (X x - 1) In ( 1 + Y^/Xj.) - (for Case 2) 


which converges for any initial guess selected from the interval 
[ 2 , X 2 ] . The suggested initial guess is 

= x 2 (for Case 2) . 


The conjugate of a given iterate is either the next or previous 
iterate, i.e., the solution is bracketed between any pair of 
adjacent iterates. 

Case 3 is defined by 


0 < X^ < 1 and 

(l-X^lnC (1-X 1 )/X 1 ] < X 2 < 2 +(1-X 1 )ln(l+X 2 /X 1 ) (defines Case 3) 

and is encountered when there is an HRR and we want to use (C2) 
to solve for P at a point close to and on either side of the ARB 
(a transitional region) . The iteration is 


Y (i+1) = Y (i) 

- (1/2) [Y^ + (1 - x x ) In (1 + Y^J/X-l) - x 2 ] (for Case 3) 


which converges for any initial guess selected from the interval 
[1/2 -X lr X 2 ]. Note that 1/2 - X ^ can be the initial guess even 
when negative, but it is not a very good guess when negative. The 
suggested initial guess is 
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(for Case 3) . 


Y (0) = X 2 

The solution is bracketed by any iterate Y (l) and its conjugate 
Y c ^) defined by 

Y c (i) s X 2 - (1 - X x ) ln(l + Y^/Xi) (for Case 3) . 

Note that the iteration can be written more concisely as 
y(i+l) = (1/2) [Y^ + Y c ^^] (for Case 3) . 

Case 4 is defined by 
0 < X x < 1 and 

0 < X 2 < (1 - X x ) In [ ( 1 - X^/X^ (defines Case 4) . 

This case will be encountered when there is a wide HRR and we 
want to use (C2) to solve for P at some point in the HRR not too 
close to the ARB. The iteration is 

Y (i+1) = (i/2) [Y< lJ - X x ] 

+ ( Xl /2) exp[(X 2 - Y^h/d - X x )] (for Case 4) 

which converges for any initial guess selected from the interval 
[0,1]. The suggested initial guess is 

y(°) = o (for Case 4) . 

The solution is bracketed by any iterate Y< 1} and its conjugate 
Y c (^) defined by 
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Y C (1) s X 1 ex P[( X 2 " Y (i) ) / (1 - X x )] - X ± 


(for Case 4) . 


Note that the iteration can be written more concisely as 
Y (i+1) = (1/2) [Yt 1 ) + Y^ 1 )] (for Case 4) . 

The following function subprogram can be appended to a FORTRAN 
source code, allowing the code to call the function F as it would 
call any built-in function. The iterations are terminated when 
error estimates indicate that the sum F(X 1 ,X 2 )+X 1 has an error 
less than one part per ten thousand. The number of iterations 
needed to produce this accuracy depends on the individual case. 
The number can be as large as twelve or thirteen (comparable to 
the bisection method) or as small as two or three. 


FUNCTION F (XI , X2 ) 

C This function subprogram can be appended to a FORTRAN source 
C code, allowing the code to call the function F defined in the 
C text. XI must be positive and X2 must be nonnegative. 

C 

DELT0L=1 . OE-4 
C 

C Check for illegal arguments. 

C 

IF (XI . LE .0.0) THEN 

WRITE (*,*) 'ERROR: XI IS NOT POSITIVE* 

GO TO 100 
END IF 

IF (X2.LT.0.0) THEN 

WRITE (*,*) 'ERROR: X2 IS NEGATIVE' 

GO TO 100 
END IF 
C 

C Determine which of the four cases apply and go to the 
C appropriate block. 

C 
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c 

c 

c 


c 

c 

c 


IF (XI . GE .1.0) GO TO 10 

XH=2 . 0+ ( 1 . 0-X1) *ALOG ( 1 . 0+X2 /XI ) 

XL= ( 1 . 0-X1) *ALOG ( (1.0-X1) /XI) 

IF (X2.GT.XH) GO TO 30 

IF ( (X2.LE.XH) .AND. (X2.GT.XL) ) GOTO 

GO TO 70 

Case 1 block starts here. 

10 CONTINUE 
Y=X2 

20 CONTINUE 

Y=X2+ (Xl-1 . 0) *ALOG ( 1 . 0+Y/X1) 

T=X2+ (Xl-1 . 0) *ALOG ( 1 . 0+Y/X1) -Y 
YC=Y+ (X2+X1) *T/ (X2+1 .0) 

ERROR=ABS ( Y-YC) / (Y+Xl) 

IF (ERROR. GT.DELTOL) GO TO 20 
GO TO 90 

Case 2 block starts here. 


50 


C 

C 

C 


30 CONTINUE 
Y=X2 

40 CONTINUE 
YC=Y 

T=X2+ (Xl-1 • 0 ) *ALOG (1 . 0+Y/X1) -Y 
Y=Y+ (X2+X1) *T / (X2+1.0) 
ERROR=ABS (Y-YC) / (Y+Xl) 

IF (ERROR. GT.DELTOL) GO TO 40 
GO TO 90 

Case 3 block starts here. 


C 

C 


50 CONTINUE 
Y=X2 

YC=X2- ( 1 . 0-X1) *ALOG ( 1 . 0+Y/X1) 
60 CONTINUE 

Y=0 . 5* ( Y+YC) 

YC=X2- ( 1 • 0-X1) *ALOG ( 1 . 0+Y/X1) 
ERROR=ABS (Y-YC) / (Y+Xl) 

IF (ERROR. GT.DELTOL) GO TO 60 
GO TO 90 

Case 4 block starts here. 
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c 

70 CONTINUE 
Y=0 . 0 

T= (X2-Y) / (1.0-X1) 

YC=X1*EXP (T) -XI 
80 CONTINUE 

Y=0 . 5* ( Y+YC) 

T= (X2-Y) / (1.0-X1) 

YC=X1*EXP (T) -XI 
ERROR=ABS (Y-YC) / ( Y+Xl) 

IF (ERROR. GT.DELTOL) GO TO 80 
90 CONTINUE 
F=Y 

100 CONTINUE 
RETURN 
END 
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